 
      subroutine lbfgsb(n,m,x,l,u,nbd,f,g,factr,wa,iwa,task,iprint,
     +                              isbmin,csave,lsave,isave,dsave)
 
      character*60 task,csave
      logical lsave(4)
      integer n,m,iprint,isbmin,
     +        nbd(n),iwa(3*n),isave(44)
      double precision f,factr,
     +       x(n),l(n),u(n),g(n),wa(2*m*n+4*n+12*m*m+12*m),dsave(29)
 
c     ************
c
c     Subroutine lbfgsb
c
c     This subroutine partitions the working arrays wa and iwa, and 
c       then uses the limited memmory BFGS method to solve the bound
c       constrained optimization problem by calling mainlb.
c
c     The subroutine statement is:
c
c       subroutine lbfgsb(n,m,x,l,u,nbd,f,g,factr,wa,iwa,task,iprint,
c                                     isbmin,csave,lsave,isave,dsave)
c
c       where
c
c     n is an integer variable.
c       On entry n is the dimension of the problem.
c       On exit n is unchanged.
c
c     m is an integer variable.
c       On entry m is the maximum number of variable metric corrections
c         used to define the limited memory matrix.
c       On exit m is unchanged.
c
c     x is a double precision array of dimension n.
c       On entry x is an approximation to the solution.
c       On exit x is the current approximation.
c
c     l is a double precision array of dimension n.
c       On entry l is the lower bound of x.
c       On exit l is unchanged
c
c     u is a double precision array of dimension n.
c       On entry u is the upper bound of x.
c       On exit u is unchanged.
c
c     nbd is an integer array of dimension n.
c       On entry nbd represents the type of bounds imposed on the
c         variables, and must be specified as follows:
c         nbd(i)=0 if x(i) is unbounded,
c                1 if x(i) has only a lower bound,
c                2 if x(i) has both lower and upper bounds, and
c                3 if x(i) has only an upper bound.
c       On exit nbd is unchanged.
c
c     f is a double precision variable.
c       On first entry f is unspecified.
c       On final exit f is the value of the function at x.
c
c     g is a double precision array of dimension n.
c       On first entry g is unspecified.
c       On final exit g is the value of the gradient at x.
c
c     factr is a double precision variable.
c       On entry factr >= 0 is specified by the user.  The iteration
c         will stop when
c
c         (f^k - fb^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch
c
c         where epsmch is the machine precision which is automatically
c         generated by the code. 
c       On exit factr is unchanged.
c
c     wa is a double presivion working array of length 
c       (2mmax + 4)nmax + 12mmax^2 + 12mmax.
c
c     iwa is an integer working array of length 3nmax.
c
c     task is a working string of characters of length 60 indicating
c       the current job when entering and quiting this subroutine.
c
c     iprint is an integer variable that must be set by the user.
c       It controls the frequency and type of output generated:
c        iprint<0    no output is generated;
c        iprint=0    print only one line at the last iteration;
c        0<iprint<99 print also f and |proj g| every iprint iterations;
c        iprint=99   print details of every iteration except n-vectors;
c        iprint=100  print also the changes of active set and final x;
c        iprint>100  print details of every iteration including x and g;
c       When iprint > 0, the file iterate.dat will be created to
c                        summarize the iteration.
c
c     isbmin is an integer variable.
c       On entry isbmin specifies the subspace minimization methods.
c         isbmin = 1  the direct method will be used,
c                  2  the dual method will be used,
c                  3  the conjugate gradient method will be used.
c       On exit isbmin is unchanged.
c
c     csave is a working string of characters of length 60.
c
c     lsave is a logical working array of dimension 4.
c       On exit with 'task' = NEW_X, the following information are
c                                                             available:
c         If lsave(1) = .true.  then  the initial X is replaced by its
c                                     projection in the feasible set;
c         If lsave(2) = .true.  then  the problem is constrained;
c         If lsave(3) = .true.  then  each variable has upper and lower
c                                     bounds;
c
c     isave is an integer working array of dimension 44.
c       On exit with 'task' = NEW_X, the following information are
c                                                             available:
c         isave(22) = the total number of intervals explored in the 
c                         search of Cauchy points;
c         isave(23) = the total number of CG iterations;
c         isave(26) = the total number of skipped BFGS update before 
c                         the current iteration;
c         isave(30) = the number of current iteration;
c         isave(31) = the total number of BFGS update prior the current
c                         iteration;
c         isave(32) = the number of CG iteration in the current
c                         iteration;
c         isave(33) = the number of intervals explored in the search of
c                         Cauchy point in the current iteration;
c         isave(34) = the total number of function and gradient 
c                         evaluations;
c         isave(36) = the number of function value or gradient
c                                  evaluations in the current iteration;
c         if isave(37) = 0  then the subspace argmin is within the box;
c         if isave(37) = 1  then the subspace argmin is beyond the box;
c         isave(38) = the number of free variables in the current
c                         iteration;
c         isave(39) = the number of active constraints in the current
c                         iteration;
c         n + 1 - isave(40) = the number of variables leaving the set of
c                            active constraint in the current iteration;
c         isave(41) = the number of variables entering the set of active
c                         constraint in the current iteration.
c
c     dsave is a double precision working array of dimension 29.
c       On exit with 'task' = NEW_X, the following information are
c                                                             available:
c         dsave(1) = current 'theta' in the BFGS matrix;
c         dsave(2) = f(x) in the previous iteration;
c         dsave(3) = factr*epsmch;
c         dsave(4) = 2-norm of the line search direction vector;
c         dsave(5) = the machine precision epsmch generated by the code;
c         dsave(7) = the accumulated time spent on searching for
c                                                         Cauchy points;
c         dsave(8) = the accumulated time spent on
c                                                 subspace minimization;
c         dsave(9) = the accumulated time spent on line search;
c         dsave(11) = the slope of the line search function at
c                                  the current point of line search;
c         dsave(12) = the maximum relative step length imposed in
c                                                           line search;
c         dsave(13) = the infinity norm of the projected gradient;
c         dsave(14) = the relative step length in line search;
c         dsave(15) = the slope of the line search function at
c                                  the starting point of line search;
c         dsave(16) = the square of the 2-norm of the line search
c                                                      direction vector.
c
c     Subprograms called:
c
c       L-BFGS-B Library ... mainlb.    
c
c
c     References:
c       R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
c       memory algorithm for bound constrained optimization'',
c       SIAM J. Scientific Computing 16 (1995), no. 5.
c
c       C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a
c       limited memory FORTRAN code for solving bound constrained
c       optimization problems'', Tech. Report, EECS Department,
c       Northwestern University, 1994.
c
c       (Postscript files of these papers are available via anonymous
c        ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************
 
      integer l1,l2,l3,lws,lr,lz,lt,ld,lsg,lwa,lyg,
     +        lsgo,lwy,lsy,lss,lyy,lwt,lwn,lsnd,lygo

      if (task .eq. 'START') then
         isave(1) = m*n
         isave(2) = m**2
         isave(3) = 4*m**2
         isave(4) = 1
         isave(5) = isave(4) + isave(1)
         isave(6) = isave(5) + isave(1)
         isave(7) = isave(6) + isave(2)
         isave(8) = isave(7) + isave(2)
         isave(9) = isave(8) + isave(2)
         isave(10) = isave(9) + isave(2)
         isave(11) = isave(10) + isave(3)
         isave(12) = isave(11) + isave(3)
         isave(13) = isave(12) + n
         isave(14) = isave(13) + n
         isave(15) = isave(14) + n
         isave(16) = isave(15) + n
         isave(17) = isave(16) + 8*m
         isave(18) = isave(17) + m
         isave(19) = isave(18) + m
         isave(20) = isave(19) + m   
      endif
         l1 = isave(1)
         l2 = isave(2)
         l3 = isave(3)
         lws = isave(4)
         lwy = isave(5)
         lsy = isave(6)
         lss = isave(7)
         lyy = isave(8)
         lwt = isave(9)
         lwn = isave(10)
         lsnd = isave(11)
         lz = isave(12)
         lr = isave(13)
         ld = isave(14)
         lt = isave(15)
         lwa = isave(16)
         lsg = isave(17)
         lsgo = isave(18)
         lyg = isave(19)
         lygo = isave(20)
      call mainlb(n,m,x,l,u,nbd,f,g,factr,
     +  wa(lws),wa(lwy),wa(lsy),wa(lss),wa(lyy),wa(lwt),
     +  wa(lwn),wa(lsnd),wa(lz),wa(lr),wa(ld),wa(lt),
     +  wa(lwa),wa(lsg),wa(lsgo),wa(lyg),wa(lygo),
     +  iwa(1),iwa(n+1),iwa(2*n+1),task,iprint,isbmin,
     +  csave,lsave,isave(22),dsave)

      return

      end

c======================= The end of lbfgsb =============================

      subroutine mainlb(n,m,x,l,u,nbd,f,g,factr,ws,wy,sy,ss,yy,wt,
     +                   wn,snd,z,r,d,t,wa,sg,sgo,yg,ygo,index,
     +                   iwhere,indx2,task,iprint,isbmin,
     +                   csave,lsave,isave,dsave)
 
      character*60 task,csave
      logical lsave(4)
      integer n,m,iprint,isbmin,
     +        nbd(n),index(n),iwhere(n),indx2(n),isave(23)
      double precision f,factr,
     +       x(n),l(n),u(n),g(n),z(n),r(n),d(n),t(n),
     +       wa(8*m),sg(m),sgo(m),yg(m),ygo(m),
     +       ws(n,m),wy(n,m),sy(m,m),ss(m,m),yy(m,m),wt(m,m),
     +       wn(2*m,2*m),snd(2*m,2*m),dsave(29)

c     ************
c
c     Subroutine mainlb
c
c     This subroutine solves bound constrained optimization problems by
c       using the compact formula of the limited memory BFGS updates.
c       
c     The subroutine statement is:
c
c     subroutine mainlb(n,m,x,l,u,nbd,f,g,factr,ws,wy,sy,ss,yy,wt,
c                        wn,snd,z,r,d,t,wa,sg,sgo,yg,ygo,index,
c                        iwhere,indx2,task,iprint,isbmin,
c                        csave,lsave,isave,dsave)
c
c       where
c
c     n is an integer variable.
c       On entry n is total number of variables.
c       On exit n is unchanged.
c
c     m is an integer variable.
c       On entry m is the maximum number of variable metric
c          corrections allowed in the limited memory matrix.
c       On exit m is unchanged.
c
c     x is a double precision array of dimension n.
c       On entry x is an approximation to the solution.
c       On exit x is the current approximation.
c
c     l is a double precision array of dimension n.
c       On entry l is the lower bound of x.
c       On exit l is unchanged.
c
c     u is a double precision array of dimension n.
c       On entry u is the upper bound of x.
c       On exit u is unchanged.
c
c     nbd is an integer array of dimension n.
c       On entry nbd represents the type of bounds imposed on the
c         variables, and must be specified as follows:
c         nbd(i)=0 if x(i) is unbounded,
c                1 if x(i) has only a lower bound,
c                2 if x(i) has both lower and upper bounds,
c                3 if x(i) has only upper bound.
c       On exit nbd is unchanged.
c
c     f is a double precision variable.
c       On first entry f is unspecified.
c       On final exit f is the value of the function at x.
c
c     g is a double precision array of dimension n.
c       On first entry g is unspecified.
c       On final exit g is the value of the gradient at x.
c
c     factr is a double precision variable. 
c       on entry factr > 1 is specified by the user.  The iteration will
c         be stopped once the relative reduction of the objective value 
c         becomes less than factr*epsmch, where epsmch is the
c         machine precision automatically generated by the code.
c       On exit factr is unchanged.
c
c     ws, wy, sy, and wt are double precision working arrays used to
c       store the following information defining the limited memory
c          BFGS matrix:
c          ws, of dimension n x m, stores S, the matrix of s-vectors;
c          wy, of dimension n x m, stores Y, the matrix of y-vectors;
c          sy, of dimension m x m, stores S'Y;
c          ss, of dimension m x m, stores S'S;
c	   yy, of dimension m x m, stores Y'Y;
c          wt, of dimension m x m, stores the Cholesky factorization
c                                  of (theta*S'S+LD^(-1)L'); see eq.
c                                  (2.26) in [1].
c
c     wn is a double precision working array of dimension 2m x 2m
c       used to store the LDL^T factorization of the indefinite matrix
c                 N = [Y' ZZ'Y   L_a'+R_z']
c                     [L_a'+R_z' S'AA'S   ]
c
c     snd is a double precision working array of dimension 2m x 2m
c       used to store the lower triangular part of
c                 N = [Y' ZZ'Y   L_a'+R_z']
c                     [L_a'+R_z' S'AA'S   ]
c	     
c     z(n),r(n),d(n),t(n),wa(8*m) are double precision working arrays.
c       z is used at different times to store the Cauchy point and
c       the Newton point.
c
c     sg(m),sgo(m),yg(m),ygo(m) are double precision working arrays. 
c
c     index is an integer working array of dimension n.
c       In subroutine freev, index is used to store the free and fixed
c          variables at the GCP.
c
c     iwhere is an integer working array of dimension n used to record
c       the status of the vector x for GCP computation.
c       iwhere(i)=0 or -3 if x(i) is free and has bounds,
c                 1       if x(i) is fixed at l(i), and l(i) .ne. u(i)
c                 2       if x(i) is fixed at u(i), and u(i) .ne. l(i)
c                 3       if x(i) is always fixed, i.e.,  u(i)=x(i)=l(i)
c                 -1      if x(i) is always free, i.e., no bounds on it.
c
c     indx2 is an integer working array of dimension n.
c       Within subroutine cauchy, indx2 corresponds to the array iorder.
c       In subroutine freev, a list of variables entering and leaving
c       the free set is stored in indx2, and it is passed on to
c       subroutine formk with this information.
c
c     task is a working string of characters of length 60 indicating
c       the current job when entering and quiting this subroutine.
c
c     iprint is an INTEGER variable that must be set by the user.
c       It controls the frequency and type of output generated:
c        iprint<0    no output is generated;
c        iprint=0    print only one line at the last iteration;
c        0<iprint<99 print also f and |proj g| every iprint iterations;
c        iprint=99   print details of every iteration except n-vectors;
c        iprint=100  print also the changes of active set and final x;
c        iprint>100  print details of every iteration including x and g;
c       When iprint > 0, the file iterate.dat will be created to
c                        summarize the iteration.
c
c     isbmin is an integer variable.
c       On entry isbmin specifies the subspace minimization methods.
c         isbmin = 1  the direct method will be used,
c                  2  the dual method will be used,
c                  3  the conjugate gradient method will be used.
c       On exit isbmin is unchanged.
c
c     csave is a working string of characters of length 60.
c
c     lsave is a logical working array of dimension 4.
c
c     isave is an integer working array of dimension 23.
c
c     dsave is a double precision working array of dimension 29.
c
c
c     Subprograms called
c
c       L-BFGS-B Library ... cauchy, subsm, subdl, subcg, lnsrch, formk, 
c
c        errchk, print1, print2, print3, active, projgr,
c
c        freev, compr, matupd, formt, sytimg, formh, heapst.
c
c       Minpack2 Library ... timer, epmeps.
c
c       Linpack Library ... dcopy, ddot.
c
c
c     References:
c       [1] R. Byrd, J. Nocedal and R. Schnabel "Representations of
c       Quasi-Newton Matrices and their use in Limited Memory Methods,
c       (1992), Northwestern University,  EECS Dept., Rep. NAM 06, to
c       appear in Mathematical Programming.
c
c       R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
c       memory algorithm for bound constrained optimization'',
c       SIAM J. Scientific Computing 16 (1995), no. 5.
c
c       C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a
c       limited memory FORTRAN code for solving bound constrained
c       optimization problems'', Tech. Report, EECS Department,
c       Northwestern University, 1994.
c
c       (Postscript files of these papers are available via anonymous
c        ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************
 
      double precision one,zero
      parameter(one=1.0d0,zero=0.0d0)

      logical prjctd,cnstnd,boxed,updatd,wrk
      character*3 word
      integer i,k,nintol,ncgtol,itfile,iback,nskip,
     +        head,col,iter,itail,iupdat,
     +        itercg,nint,nfgv,info,ifun,
     +        iword,nfree,nact,ileave,nenter
      double precision theta,fold,ddot,dr,rr,tol,dpmeps,
     +                 xstep,sbgnrm,ddum,dnorm,dtd,epsmch,
     +                 cpu1,cpu2,cachyt,sbtime,lnscht,time1,time2,
     +                 gd,gdold,stp,stpmx,time
      
      if (task .eq. 'START') then

         call timer(time1)

c        Generate the current machine precision.

         epsmch = dpmeps()

c        Initialize counters and scalars when task='START'.

c           for the limited memory BFGS matrices:
         col = 0
         head = 1
         theta = one
         iupdat = 0
         updatd = .false.
 
c           for operation count:
         iter = 0
         nfgv = 0
         nint = 0
         nintol = 0
         ncgtol = 0
         nskip = 0
         nfree = n

c           tolerance in stop criterion.
         tol = factr*epsmch

c           for measuring running time:
         cachyt = 0
         sbtime = 0
         lnscht = 0
 
c           'word' records the status of subspace solutions.
         word = '---'

c           'info' records the termination information.
         info = 0

         if (iprint .ge. 1) then
c                                open a summary file 'iterate.dat'
            open (8, file = 'iterate.dat', status = 'unknown')
            itfile = 8
         endif            

c        Check the input arguments for errors.

	 call errchk(n,m,factr,isbmin,l,u,nbd,task,info,k)
         if (task(1:5) .eq. 'ERROR') then
            call print3(n,l,u,nbd,x,f,task,iprint,info,itfile,
     +                  iter,nfgv,nintol,ncgtol,nskip,nact,sbgnrm,
     +                  zero,nint,word,iback,stp,xstep,k,
     +                  cachyt,sbtime,lnscht)
            return
         endif

         call print1(n,m,l,u,x,iprint,isbmin,itfile,epsmch)
 
c        Initialize iwhere & project the initial x to the feasible set.
 
         call active(n,l,u,nbd,x,iwhere,iprint,prjctd,cnstnd,boxed) 

c     The end of the start block.

      else
c          restore local variables.

         prjctd = lsave(1)
         cnstnd = lsave(2)
         boxed = lsave(3)
         updatd = lsave(4)

         nintol = isave(1)
         ncgtol = isave(2)
         itfile = isave(3)
         iback = isave(4)
         nskip = isave(5)
         head = isave(6)
         col = isave(7)
         itail = isave(8)
         iter = isave(9)
         iupdat = isave(10)
         itercg = isave(11)
         nint = isave(12)
         nfgv = isave(13)
         info = isave(14)
         ifun = isave(15)
         iword = isave(16)
         nfree = isave(17)
         nact = isave(18)
         ileave = isave(19)
         nenter = isave(20)

         theta = dsave(1)
         fold = dsave(2)
         tol = dsave(3)
         dnorm = dsave(4)
         epsmch = dsave(5)
         cpu1 = dsave(6)
         cachyt = dsave(7)
         sbtime = dsave(8)
         lnscht = dsave(9)
         time1 = dsave(10)
         gd = dsave(11)
         stpmx = dsave(12)
         sbgnrm = dsave(13)
         stp = dsave(14)
         gdold = dsave(15)
         dtd = dsave(16)
   
c        If f(x) and grad(x) are calculated after returning to the 
c           driver, go directly to the section where they are needed.

         if (task(1:5) .eq. 'FG_LN') goto 666
         if (task(1:5) .eq. 'NEW_X') goto 777
         if (task(1:5) .eq. 'FG_ST') goto 111
         if (task(1:4) .eq. 'STOP') then
            if (task(7:9) .eq. 'CPU') then
c                                          restore the previous iterate.
               call dcopy(n,t,1,x,1)
               call dcopy(n,r,1,g,1)
               f = fold
            endif
            goto 999
         endif
      endif 

c     Compute f0 and g0.

      task = 'FG_START' 
c          return to the driver for calculating f and g; reenter at 111.
      goto 1000
 111  continue
      nfgv = 1
 
c     Compute the infinity norm of the projected (-)gradient.
 
      call projgr(n,l,u,nbd,x,g,sbgnrm)
  
      if (iprint .ge. 1) then
         write (6,1002) iter,f,sbgnrm
         write (itfile,1003) iter,nfgv,sbgnrm,f
      endif
      if (sbgnrm .le. zero) then
c                                terminate the algorithm.
         task = 'CONVERGENCE: NORM OF PROJECTED GRADIENT = 0'
         goto 999
      endif 
 
c ----------------- the beginning of the loop --------------------------
 
 222  continue
      if (iprint .ge. 99) write (6,1001) iter + 1
      iword = -1
c
      if (.not. cnstnd .and. col .gt. 0) then 
c                                            skip the search for GCP.
         call dcopy(n,x,1,z,1)
	 wrk = updatd
         nint = 0
         goto 333
      endif

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     Compute the GCP.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      call timer(cpu1) 
      call cauchy(n,x,l,u,nbd,g,indx2,iwhere,t,d,z,
     +            m,wy,ws,sy,wt,theta,col,head,
     + 		  wa(1),wa(2*m+1),wa(4*m+1),wa(6*m+1),nint,
     +            isbmin,sg,yg,iprint,sbgnrm,info)
      if (info .ne. 0) then 
c          singular trangular system detected; refresh the lbfgs memory.
         if(iprint .ge. 1) write (6, 1005)
         info = 0
         col = 0
         head = 1
         theta = one
         iupdat = 0
         updatd = .false.
         call timer(cpu2) 
         cachyt = cachyt + cpu2 - cpu1
         goto 222
      endif
      call timer(cpu2) 
      cachyt = cachyt + cpu2 - cpu1
      nintol = nintol + nint

c     Count the entering and leaving variables for iter > 0; 
c     find the index set of free and active variables at the GCP.

      call freev(n,nfree,index,nenter,ileave,indx2,
     +           iwhere,wrk,updatd,cnstnd,iprint,iter)

      nact = n - nfree
 
 333  continue
 
c     If there is no free variables or B=I, then
c                                        skip the subspace minimization.
 
      if (nfree .eq. 0 .or. col .eq. 0) goto 555
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     Subspace minimization using the method chosen.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      call timer(cpu1) 

c     Form  the LEL^T factorization of the indefinite
c       matrix    K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
c                     [L_a -R_z           theta*S'AA'S ]
c                                                    where E = [-I  0]
c                                                              [ 0  I]

      if (isbmin .ne. 3) then
         if (wrk) call formk(n,nfree,index,nenter,ileave,indx2,iupdat,
     +                    updatd,wn,snd,m,ws,wy,sy,theta,col,head,info)
         if (info .ne. 0) then
c             nonpositive definiteness in Cholesky factorization;
c             refresh the lbfgs memory and restart the iteration.
            if(iprint .ge. 1) write (6, 1006)
            info = 0
            col = 0
            head = 1
            theta = one
            iupdat = 0
            updatd = .false.
            call timer(cpu2) 
            sbtime = sbtime + cpu2 - cpu1 
            goto 222
         endif
      endif 

      if(isbmin .eq. 2) then
c                            call the dual method.
         call subdl(n,m,nact,index,l,u,nbd,x,z,g,d,t,
     +                ws,wy,sy,yy,theta,col,head,sg,yg,iword,
     +                wa,wa(2*m+1),wn,iprint,info)
      else 
c           compute r=-Z'B(xcp-xk)-Z'g (using wa(2m+1)=W'(xcp-x)
c                                                      from 'cauchy').
         call compr(n,m,x,f,g,ws,wy,sy,wt,z,r,d,wa,index,
     +              theta,col,head,nfree,cnstnd,info)
         if (info .ne. 0) goto 444
         if(isbmin .eq. 1)then
c                              call the direct method.
            call subsm(n,m,nfree,index,l,u,nbd,z,r,ws,wy,theta,
     +                 col,head,iword,wa,wn,iprint,info)
         else
c             call the CG method.
            call subcg(n,nfree,index,l,u,nbd,z,r,d,t,m,ws,wy,sy,wt,
     +                 theta,col,head,itercg,wa,iword,iprint,info)
            ncgtol = ncgtol + itercg
	 endif
      endif
 444  continue
      if (info .ne. 0) then 
c          singular trangular system detected;
c          refresh the lbfgs memory and restart the iteration.
         if(iprint .ge. 1) write (6, 1005)
         info = 0
         col = 0
         head = 1
         theta = one
         iupdat = 0
         updatd = .false.
         call timer(cpu2) 
         sbtime = sbtime + cpu2 - cpu1 
         goto 222
      endif
 
      call timer(cpu2) 
      sbtime = sbtime + cpu2 - cpu1 
 555  continue
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     Line search and optimality tests.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
c     Generate the search direction d:=z-x.

      do 40 i = 1, n
	 d(i) = z(i) - x(i)
  40  continue
      call timer(cpu1) 
 666  continue
      call lnsrch(n,l,u,nbd,x,f,fold,gd,gdold,g,d,r,t,z,stp,dnorm,
     +            dtd,xstep,stpmx,iter,ifun,iback,nfgv,info,task,
     +            boxed,cnstnd,csave,isave(22),dsave(17))
      if (info .ne. 0 .or. iback .ge. 20) then
         if (col .eq. 0) then
c                             restore the previous iterate.
            call dcopy(n,t,1,x,1)
            call dcopy(n,r,1,g,1)
            f = fold
c             abnormal termination.
            if (info .eq. 0) then
               info = -9
c                restore the actual number of f and g evaluations etc.
               nfgv = nfgv - 1
               ifun = ifun - 1
               iback = iback - 1
            endif
            task = 'ABNORMAL_TERMINATION_IN_LNSRCH'
            iter = iter + 1
            goto 999
         else
c             refresh the lbfgs memory and restart the iteration.
            if(iprint .ge. 1) write (6, 1008)
            if (info .eq. 0) nfgv = nfgv - 1
            info = 0
            col = 0
            head = 1
            theta = one
            iupdat = 0
            updatd = .false.
            task = 'RESTART_FROM_LNSRCH'
            call timer(cpu2)
            lnscht = lnscht + cpu2 - cpu1
            goto 222
         endif
      else if (task(1:5) .eq. 'FG_LN') then
c          return to the driver for calculating f and g; reenter at 666.
	 goto 1000
      else 
c          calculate and print out the quatities related to the new X.
         call timer(cpu2) 
         lnscht = lnscht + cpu2 - cpu1
         iter = iter + 1
 
c        Compute the infinity norm of the projected (-)gradient.
 
         call projgr(n,l,u,nbd,x,g,sbgnrm)
 
c        Print iteartion information.

         call print2(n,x,f,g,iprint,itfile,iter,nfgv,nact,
     +               sbgnrm,nint,word,iword,iback,stp,xstep)
         goto 1000
      endif
 777  continue

c     Test for termination.

      if (sbgnrm .le. zero) then
c                                terminate the algorithm.
         task = 'CONVERGENCE: NORM OF PROJECTED GRADIENT = 0'
         goto 999
      endif 

      ddum = max(abs(fold), abs(f), one)
      if ((fold - f) .le. tol*ddum) then
c                                        terminate the algorithm.
         task = 'CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH'
         if (iback .ge. 10) info = -5
c           i.e., to issue a warning if iback>10 in the line search.
         goto 999
      endif 

c     Compute d=newx-oldx, r=newg-oldg, rr=y'y and dr=y's.
 
      do 42 i = 1, n
         r(i) = g(i) - r(i)
  42  continue
      rr = ddot(n,r,1,r,1)
      if (stp .eq. one) then  
         dr = gd - gdold
         ddum = -epsmch*gdold
      else
         dr = (gd - gdold)*stp
	 call dscal(n,stp,d,1)
         ddum = -epsmch*gdold*stp
      endif
 
      if (dr .le. ddum) then
c                            skip the L-BFGS update.
         nskip = nskip + 1
         updatd = .false.
         if (iprint .ge. 1) write (6,1004) dr, ddum
         goto 888
      endif 
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     Update the L-BFGS matrix.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
      updatd = .true.
      iupdat = iupdat + 1

c     Update matrices WS and WY and form the middle matrix in B.

      call matupd(n,m,ws,wy,sy,ss,d,r,itail,
     +            iupdat,col,head,theta,rr,dr,stp,dtd)

c     Form the upper half of the pds T = theta*SS + L*D^(-1)*L';
c        Store T in the upper triangular of the array wt;
c        Cholesky factorize T to J*J' with
c           J' stored in the upper triangular of wt.

      call formt(n,m,wt,sy,ss,col,theta,info)
 
      if (info .ne. 0) then 
c          nonpositive definiteness in Cholesky factorization;
c          refresh the lbfgs memory and restart the iteration.
         if(iprint .ge. 1) write (6, 1007)
         info = 0
         col = 0
         head = 1
         theta = one
         iupdat = 0
         updatd = .false.
         goto 222
      endif

c     Now the inverse of the middle matrix in B is

c       [  D^(1/2)      O ] [ -D^(1/2)  D^(-1/2)*L' ]
c       [ -L*D^(-1/2)   J ] [  0        J'          ]

 888  continue
      if (isbmin .eq. 2) then
c                             compute S'g and Y'g.
         call sytimg(n,m,g,ws,wy,sg,sgo,yg,ygo,col,head,iupdat)

         if (.not. updatd) goto 222
c                                    to skip the following update. 
 
c        Form the middle matrix in H, update YY
c                                        and the upper triangular of SY.
         call formh(n,m,sy,yy,sg,yg,sgo,ygo,rr,col,iupdat)
 
      endif
 
c -------------------- the end of the loop -----------------------------
 
      goto 222
 999  continue
      call timer(time2)
      time = time2 - time1
      call print3(n,l,u,nbd,x,f,task,iprint,info,itfile,
     +            iter,nfgv,nintol,ncgtol,nskip,nact,sbgnrm,
     +            time,nint,word,iback,stp,xstep,k,
     +            cachyt,sbtime,lnscht)
 1000 continue

c     Save local variables.

      lsave(1) = prjctd
      lsave(2) = cnstnd
      lsave(3) = boxed
      lsave(4) = updatd

      isave(1) = nintol 
      isave(2) = ncgtol 
      isave(3) = itfile 
      isave(4) = iback 
      isave(5) = nskip 
      isave(6) = head 
      isave(7) = col 
      isave(8) = itail 
      isave(9) = iter 
      isave(10) = iupdat 
      isave(11) = itercg 
      isave(12) = nint 
      isave(13) = nfgv 
      isave(14) = info 
      isave(15) = ifun 
      isave(16) = iword 
      isave(17) = nfree 
      isave(18) = nact 
      isave(19) = ileave 
      isave(20) = nenter 

      dsave(1) = theta 
      dsave(2) = fold 
      dsave(3) = tol 
      dsave(4) = dnorm 
      dsave(5) = epsmch 
      dsave(6) = cpu1 
      dsave(7) = cachyt 
      dsave(8) = sbtime 
      dsave(9) = lnscht 
      dsave(10) = time1 
      dsave(11) = gd 
      dsave(12) = stpmx 
      dsave(13) = sbgnrm
      dsave(14) = stp
      dsave(15) = gdold
      dsave(16) = dtd  

 1001 format (//,'ITERATION ',i5)
 1002 format
     +  (/,'At iterate',i5,4x,'f= ',1p,d12.5,4x,'|proj g|= ',1p,d12.5)
 1003 format (2(1x,i4),5x,'-',5x,'-',3x,'-',5x,'-',5x,'-',8x,'-',3x,
     +        1p,2(1x,d10.3))
 1004 format ('  ys=',1p,e10.3,'  -gs=',1p,e10.3,' BFGS update SKIPPED')
 1005 format (/, 
     +' Singular trangular system detected;',/,
     +'   refresh the lbfgs memory and restart the iteration.')
 1006 format (/, 
     +' Nonpositive definiteness in Cholesky factorization in formk;',/,
     +'   refresh the lbfgs memory and restart the iteration.')
 1007 format (/, 
     +' Nonpositive definiteness in Cholesky factorization in formt;',/,
     +'   refresh the lbfgs memory and restart the iteration.')
 1008 format (/, 
     +' Bad direction in the line search;',/,
     +'   refresh the lbfgs memory and restart the iteration.')

      return   

      end
 
c======================= The end of mainlb =============================

      subroutine active(n,l,u,nbd,x,iwhere,iprint,prjctd,cnstnd,boxed)

      logical prjctd,cnstnd,boxed
      integer n,iprint,nbd(n),iwhere(n)
      double precision x(n),l(n),u(n)

c     ************
c
c     Subroutine active
c
c     This subroutine initializes iwhere and projects the initial x to
c       the feasible set if necessary.
c
c     iwhere is an integer array of dimension n.
c       On entry iwhere is unspecified.
c       On exit iwhere(i)=-1  if x(i) has no bounds
c                         3   if l(i)=u(i)
c                         0   otherwise.
c       In cauchy, iwhere is given finer gradations.
c
c     nact is an integer variable. 
c       On entry nact is unspecified.
c       On exit nact is the number of variables equal to a bound.
c 
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************

      double precision zero
      parameter(zero=0.0d0)

      integer nbdd,i

c     Initialize nbdd, prjctd, cnstnd and boxed.

      nbdd = 0
      prjctd = .false.
      cnstnd = .false.
      boxed = .true.

c     Project the initial x to the easible set if necessary.

      do 10 i = 1, n
         if (nbd(i) .gt. 0) then
            if (nbd(i) .le. 2 .and. x(i) .le. l(i)) then
	       if (x(i) .lt. l(i)) then
                  prjctd = .true.
	          x(i) = l(i)
               endif
               nbdd = nbdd + 1
            else if (nbd(i) .ge. 2 .and. x(i) .ge. u(i)) then
	       if (x(i) .gt. u(i)) then
                  prjctd = .true.
	          x(i) = u(i)
               endif
               nbdd = nbdd + 1
            endif
         endif
  10  continue

c     Initialize iwhere and assign values to cnstnd and boxed.

      do 20 i = 1, n
         if (nbd(i) .ne. 2) boxed = .false.
         if (nbd(i) .eq. 0) then
c                                this variable is always free
	    iwhere(i) = -1

c           otherwise set x(i)=mid(x(i), u(i), l(i)).
         else
	    cnstnd = .true.
            if (nbd(i) .eq. 2 .and. u(i) - l(i) .le. zero) then
c                   this variable is always fixed
	       iwhere(i) = 3
            else 
               iwhere(i) = 0
            endif
         endif
  20  continue

      if (iprint .ge. 0) then
         if (prjctd) write (6,*)
     +   'The initial X is infeasible.  Restart with its projection.'
         if (.not. cnstnd)
     +      write (6,*) 'This problem is unconstrained.'
      endif

      if (iprint .gt. 0) write (6,1001) nbdd

 1001 format (/,'At X0 ',i9,' variables are exactly at the bounds') 

      return

      end

c======================= The end of active =============================
 
      subroutine bmv(m,sy,wt,col,v,p,info)

      integer m,col,info
      double precision sy(m,m),wt(m,m),v(2*col),p(2*col)

c     ************
c
c     Subroutine hmv
c
c     This subroutine computes the product of the 2m x 2m middle matrix 
c	in the compact L-BFGS formula of B and a 2m vector v;  
c	it returns the product in p.
c	
c     The subroutine statement is:
c
c       subroutine bmv(m,sy,wt,col,v,p)
c
c       where
c
c     m is an integer variable.
c       On entry m is the maximum number of variable metric corrections
c         used to define the limited memory matrix.
c       On exit m is unchanged.
c
c     sy is a double precision array of dimension m x m.
c       On entry sy specifies the matrix S'Y.
c       On exit sy is unchanged.
c
c     wt is a double precision array of dimension m x m.
c       On entry wt specifies the upper triangular matrix J' which is 
c         the Cholesky factorization of (thetaS'S+LD^(-1)L').
c       On exit wy is unchanged.
c
c     col is an integer variabl.
c       On entry col specifies the number of s-vectors (or y-vectors)
c         stored in the compact L-BFGS formula.
c       On exit col is unchanged.
c
c     v is a double precision array of dimension 2col.
c       On entry v specifies vector v.
c       On exit v is unchanged.
c
c     p is a double precision array of dimension 2col.
c       On entry p is unspecified.
c       On exit p is the product Mv.
c
c     info is an integer variable.
c       On entry info is unspecified.
c       On exit info = 0       for normal return,
c                    = nonzero for abnormal return when the system
c                                     in call of dtrsl is singular.
c
c     Subprograms called:
c
c       Linpack ... dtrsl.
c
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************
 
      integer i,k,i2
      double precision sum
 
      if (col .eq. 0) return
 
c     PART I: solve [  D^(1/2)      O ] [ p1 ] = [ v1 ]
c                   [ -L*D^(-1/2)   J ] [ p2 ]   [ v2 ].

c 	solve Jp2=v2+LD^(-1)v1.
      p(col + 1) = v(col + 1)
      do 20 i = 2, col
         i2 = col + i
         sum = 0.0d0
         do 10 k = 1, i - 1
            sum = sum + sy(i,k)*v(k)/sy(k,k)
  10     continue
         p(i2) = v(i2) + sum
  20  continue  
      call dtrsl(wt,m,col,p(col+1),11,info)
      if (info .ne. 0) return
 
c     	solve D^(1/2)p1=v1.
      do 30 i = 1, col
         p(i) = v(i)/sqrt(sy(i,i))
  30  continue 
 
c     PART II: solve [ -D^(1/2)   D^(-1/2)*L'  ] [ p1 ] = [ p1 ]
c                    [  0         J'           ] [ p2 ]   [ p2 ]. 
 
c       solve J^Tp2=p2. 
      call dtrsl(wt,m,col,p(col+1),01,info)
      if (info .ne. 0) return
 
c       compute p1=-D^(-1/2)(p1-D^(-1/2)L'p2)
c                 =-D^(-1/2)p1+D^(-1)L'p2.  
      do 40 i = 1, col
         p(i) = -p(i)/sqrt(sy(i,i))
  40  continue
      do 60 i = 1, col
         sum = 0.d0
         do 50 k = i + 1, col
            sum = sum + sy(k,i)*p(col+k)/sy(i,i)
  50     continue
         p(i) = p(i) + sum
  60  continue

      return

      end

c======================== The end of bmv ===============================

      subroutine cauchy(n,x,l,u,nbd,g,iorder,iwhere,t,d,xcp,
     +                  m,wy,ws,sy,wt,theta,col,head,p,c,wbp,
     +                  v,nint,isbmin,sg,yg,iprint,sbgnrm,info)
      
      integer n,m,head,col,nint,isbmin,iprint,info,
     +        nbd(n),iorder(n),iwhere(n)
      double precision theta,
     +       x(n),l(n),u(n),g(n),t(n),d(n),xcp(n),sg(m),yg(m),
     +       wy(n,col),ws(n,col),sy(m,m),wt(m,m),
     +       p(2*m),c(2*m),wbp(2*m),v(2*m)

c     ************
c
c     Subroutine cauchy
c
c     For given x, l, u, g (with sbgnrm > 0), and a limited memory
c       BFGS matrix B defined in terms of matrices WY, WS, WT, and
c       scalars head, col, and theta, this subroutine computes the
c       generalized Cauchy point (GCP), defined as the first local
c       minimizer of the quadratic
c
c                  Q(x + s) = g's + 1/2 s'Bs
c
c       along the projected gradient direction P(x-tg,l,u),
c       and returns the GCP in xcp. 
c       
c     The subroutine statement is:
c
c       subroutine cauchy(n,x,l,u,nbd,g,iorder,iwhere,t,d,xcp,
c                         m,wy,ws,sy,wt,theta,col,head,p,c,wbp,
c                         v,nint,isbmin,sg,yg,iprint,sbgnrm,info)
c
c       where
c
c     n is an integer variable.
c       On entry n is the dimension of the problem.
c       On exit n is unchanged.
c
c     x is a double precision array of dimension n.
c       On entry x is the starting point for the GCP computation.
c       On exit x is unchanged.
c
c     l is a double precision array of dimension n.
c       On entry l is the lower bound of x.
c       On exit l is unchanged.
c
c     u is a double precision array of dimension n.
c       On entry u is the upper bound of x.
c       On exit u is unchanged.
c
c     nbd is an integer array of dimension n.
c       On entry nbd represents the type of bounds imposed on the
c         variables, and must be specified as follows:
c         nbd(i)=0 if x(i) is unbounded,
c                1 if x(i) has only lower bound,
c                2 if x(i) has both lower and upper bounds, and
c                3 if x(i) has only upper bound. 
c       On exit nbd is unchanged.
c
c     g is a double precision array of dimension n.
c       On entry g is the gradient of f(x).  g must be a nonzero vector.
c       On exit g is unchanged.
c
c     iorder is an integer working array of dimension n.
c       iorder will be used to store the breakpoints in the piecewise
c       linear path and free variables encountered. On exit,
c         iorder(1),...,iorder(nleft) are indices of breakpoints
c                                which have not been encountered; 
c         iorder(nleft+1),...,iorder(nbreak) are indices of
c                                     encountered breakpoints; and
c         iorder(nfree),...,iorder(n) are indices of variables which
c                 have no bound constraits along the search direction.
c
c     iwhere is an integer array of dimension n.
c       On entry iwhere indicates only the permanently fixed (iwhere=3)
c       or free (iwhere= -1) components of x.
c       On exit iwhere records the status of the current x variables.
c       iwhere(i)=-3  if x(i) is free and has bounds, but is not moved
c                 0   if x(i) is free and has bounds, and is moved
c                 1   if x(i) is fixed at l(i), and l(i) .ne. u(i)
c                 2   if x(i) is fixed at u(i), and u(i) .ne. l(i)
c                 3   if x(i) is always fixed, i.e.,  u(i)=x(i)=l(i)
c                 -1  if x(i) is always free, i.e., no bounds on it.
c
c     t is a double precision working array of dimension n. 
c       t will be used to store the break points.
c
c     d is a double precision working array of dimension n.
c       d will be used to store the Cauchy direction P(x-tg)-x.
c
c     xcp is a double precision working array of dimension n.
c       xcp will be used to return the GCP on exit.
c
c     m is an integer variable.
c       On entry m is the maximum number of variable metric corrections 
c         used to define the limited memory matrix.
c       On exit m is unchanged.
c
c     ws, wy, sy, and wt are double precision arrays.
c       On entry they store information that defines the
c                             limited memory BFGS matrix:
c         ws(n,m) stores S, a set of s-vectors;
c         wy(n,m) stores Y, a set of y-vectors;
c         sy(m,m) stores S'Y;
c         wt(m,m) stores the
c                 cholesky factorization of (theta*S'S+LD^(-1)L').
c       On exit they are unchanged.
c
c     theta is a double precision variable.
c       On entry theta is the scaling factor specifying B_0 = theta I.
c       On exit theta is unchanged.
c
c     col is an integer variable.
c       On entry col is the actual number of variable metric
c         corrections stored so far.
c       On exit col is unchanged.
c
c     head is an integer variable.
c       On entry head is the location of the first s-vector (or y-vector)
c         in S (or Y).
c       On exit col is unchanged.
c
c     p is a double precision working array of dimension 2m.
c       p will be used to store the vector p = W^(T)d.
c
c     c is a double precision working array of dimension 2m.
c       c will be used to store the vector c = W^(T)(xcp-x).
c
c     wbp is a double precision working array of dimension 2m.
c       wbp will be used to store the row of W corresponding
c         to a breakpoint.
c
c     v is a double precision working array of dimension 2m.
c
c     nint is an integer variable.
c       On exit nint records the number of quadratic segments explored
c         in searching for the GCP.
c
c     isbmin is an integer variable.
c       On entry isbmin specifies the subspace minimization methods.
c         isbmin = 1  the direct method will be used,
c                  2  the dual method will be used,
c                  3  the conjugate gradient method will be used.
c       On exit isbmin is unchanged.
c
c     sg and yg are double precision arrays of dimension m.
c       On entry sg  and yg store S'g and Y'g correspondingly.
c       On exit they are unchanged. 
c 
c     iprint is an INTEGER variable that must be set by the user.
c       It controls the frequency and type of output generated:
c        iprint<0    no output is generated;
c        iprint=0    print only one line at the last iteration;
c        0<iprint<99 print also f and |proj g| every iprint iterations;
c        iprint=99   print details of every iteration except n-vectors;
c        iprint=100  print also the changes of active set and final x;
c        iprint>100  print details of every iteration including x and g;
c       When iprint > 0, the file iterate.dat will be created to
c                        summarize the iteration.
c
c     sbgnrm is a double precision variable.
c       On entry sbgnrm is the norm of the projected gradient at x.
c       On exit sbgnrm is unchanged.
c
c     info is an integer variable.
c       On entry info is 0.
c       On exit info = 0       for normal return,
c                    = nonzero for abnormal return when the system
c                                      in call of bmv is singular.
c
c     Subprograms called:
c 
c       L-BFGS-B Library ... heapst, bmv.
c
c       Linpack ... dscal dcopy, daxpy.
c
c
c     References:
c       R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
c       memory algorithm for bound constrained optimization'',
c       SIAM J. Scientific Computing 16 (1995), no. 5.
c
c       C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a
c       limited memory FORTRAN code for solving bound constrained
c       optimization problems'', Tech. Report, EECS Department,
c       Northwestern University, 1994.
c
c       (Postscript files of these papers are available via anonymous
c        ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************

      double precision one,zero
      parameter(one=1.0d0,zero=0.0d0)

      logical xlower,xupper,bnded
      integer i,j,col2,nfree,nbreak,pointr,ibp,nleft,ibkmin,iter
      double precision f1,f2,dt,dtm,tsum,dibp,zibp,dibp2,bkmin,
     +                 tu,tl,wmc,wmp,wmw,ddot,tj,tj0,neggi,sbgnrm
 
c     Check the status of the variables, reset iwhere(i) if necessary;
c       compute the Cauchy direction d and the breakpoints t; initialize
c       the derivative f1 and the vector p = W'd (for theta = 1).
 
      if (sbgnrm .le. zero) then
         if (iprint .ge. 0) write (6,*) 'Subgnorm = 0.  GCP = X.'
         call dcopy(n,x,1,xcp,1)
	 return
      endif 
      bnded = .true.
      nfree = n + 1
      nbreak = 0
      ibkmin = 0
      bkmin = zero
      col2 = 2*col
      f1 = zero
      if (iprint .ge. 99) write (6,3010)

c     If dual method is used we have -W'g available to initialize p.
c       Otherwise, we set p to zero and build it up as we determine d.

      if (isbmin .eq. 2) then
         do 10 i = 1, col
            p(i) = -yg(i)
            p(col + i) = -sg(i)
  10     continue 
      else
         do 20 i = 1, col2
            p(i) = zero
  20     continue 
      endif 

c     In the following loop we determine for each variable its bound
c        status, and its breakpoint, and update p accordingly.
c        Smallest breakpoint is identified.

      do 50 i = 1, n 
         neggi = -g(i)      
         if (iwhere(i) .ne. 3 .and. iwhere(i) .ne. -1) then
c             if x(i) is not a constant and has bounds,
c             compute the difference between x(i) and its bounds.
            if (nbd(i) .le. 2) tl = x(i) - l(i)
            if (nbd(i) .ge. 2) tu = u(i) - x(i)

c           If a variable is close enough to a bound
c             we treat it as at bound.
            xlower = nbd(i) .le. 2 .and. tl .le. zero
            xupper = nbd(i) .ge. 2 .and. tu .le. zero

c              reset iwhere(i).
            iwhere(i) = 0
            if (xlower) then
               if (neggi .le. zero) iwhere(i) = 1
            else if (xupper) then
               if (neggi .ge. zero) iwhere(i) = 2
            else
               if (abs(neggi) .le. zero) iwhere(i) = -3
            endif
         endif 
         pointr = head
         if (iwhere(i) .ne. 0 .and. iwhere(i) .ne. -1) then
            d(i) = zero
            if (isbmin .eq. 2 .and. neggi .ne. zero) then
c                                                p := p + W'e_i* (g_i).
               do 30 j = 1, col
                  p(j) = p(j) -  wy(i,pointr)* neggi
                  p(col + j) = p(col + j) - ws(i,pointr)*neggi
                  pointr = mod(pointr,m) + 1
  30           continue 
            endif
         else
            d(i) = neggi
            f1 = f1 - neggi*neggi
            if (isbmin .ne. 2) then
c                                   p := p - W'e_i* (g_i).
               do 40 j = 1, col
                  p(j) = p(j) +  wy(i,pointr)* neggi
                  p(col + j) = p(col + j) + ws(i,pointr)*neggi
                  pointr = mod(pointr,m) + 1
  40           continue 
	    endif
            if (nbd(i) .le. 2 .and. nbd(i) .ne. 0
     +                        .and. neggi .lt. zero) then
c                                 x(i) + d(i) is bounded; compute t(i).
               nbreak = nbreak + 1
               iorder(nbreak) = i
               t(nbreak) = tl/(-neggi)
	       if (nbreak .eq. 1 .or. t(nbreak) .lt. bkmin) then
		  bkmin = t(nbreak)
		  ibkmin = nbreak
               endif
            else if (nbd(i) .ge. 2 .and. neggi .gt. zero) then
c                                 x(i) + d(i) is bounded; compute t(i).
               nbreak = nbreak + 1
               iorder(nbreak) = i
               t(nbreak) = tu/neggi
	       if (nbreak .eq. 1 .or. t(nbreak) .lt. bkmin) then
		  bkmin = t(nbreak)
		  ibkmin = nbreak
               endif
            else
c                x(i) + d(i) is not bounded.
               nfree = nfree - 1
               iorder(nfree) = i
               if (abs(neggi) .gt. zero) bnded = .false.
            endif
         endif
  50  continue 
 
c     Indices of the nonzeros components of the vector d are now stored
c       in iorder(1),...,iorder(nbreak) and iorder(nfree),...,iorder(n).
c       The smallest of the nbreak breakpoints is in t(ibkmin)=bkmin.
 
      if (theta .ne. one) then
c                   complete the initialization of p for theta not= one.
         call dscal(col,theta,p(col+1),1)
      endif
 
c     Initialize GCP xcp = x.

      call dcopy(n,x,1,xcp,1)

      if (nbreak .eq. 0 .and. nfree .eq. n + 1) then
c                  is a zero vector, return with the initial xcp as GCP.
         if (iprint .gt. 100) write (6,1010) (xcp(i), i = 1, n)
         return
      endif    
 
c     Initialize c = W'(xcp - x) = 0.
  
      do 60 j = 1, col2
         c(j) = zero
  60  continue 
 
c     Initialize derivative f2.
 
      f2 =  -theta*f1 
      if (col .gt. 0) then
     	 call bmv(m,sy,wt,col,p,v,info)
	 if (info .ne. 0) return
     	 f2 = f2 - ddot(col2,v,1,p,1)
      endif
      dtm = -f1/f2
      tsum = zero
      nint = 1
      if (iprint .ge. 99) 
     +   write (6,*) 'There are ',nbreak,'  breakpoints '
 
c     If there are no breakpoints, locate the GCP and return. 
 
      if (nbreak .eq. 0) goto 888
             
      nleft = nbreak
      iter = 1
 
 
      tj = zero
 
c------------------- the beginning of the loop -------------------------
 
 777  continue
 
c     Find the next smallest breakpoint;
c       compute dt = t(nleft) - t(nleft + 1).
 
      tj0 = tj
      if (iter .eq. 1) then
c         Since we already have the smallest breakpoint we need not do
c         heapsort yet. Often only one breakpoint is used and the
c         cost of heapsort is avoided.
	 tj = bkmin
	 ibp = iorder(ibkmin)
      else
         if (iter .eq. 2) then
c             Replace the already used smallest breakpoint with the
c             breakpoint numbered nbreak > nlast, before heapsort call.
            if (ibkmin .ne. nbreak) then
               t(ibkmin) = t(nbreak)
	       iorder(ibkmin) = iorder(nbreak)
            endif 
c        Update heap structure of breakpoints
c           (if iter=2, initialize heap).
         endif
         call heapst(nleft,t,iorder,iter-2)
         tj = t(nleft)
         ibp = iorder(nleft)  
      endif 
	 
      dt = tj - tj0
 
      if (dt .ne. zero .and. iprint .ge. 100) then
         write (6,4011) nint,f1,f2
         write (6,5010) dt
         write (6,6010) dtm
      endif	     
 
c     If a minimizer is within this interval, locate the GCP and return. 
 
      if (dtm .lt. dt) goto 888
 
c     Otherwise fix one variable and
c       reset the corresponding component of d to zero.
    
      tsum = tsum + dt
      nleft = nleft - 1
      iter = iter + 1
      dibp = d(ibp)
      d(ibp) = zero
      if (dibp .gt. zero) then
	 zibp = u(ibp) - x(ibp)
	 xcp(ibp) = u(ibp)
         iwhere(ibp) = 2
      else
	 zibp = l(ibp) - x(ibp)
	 xcp(ibp) = l(ibp)
         iwhere(ibp) = 1
      endif
      if (iprint .ge. 100) write (6,*) 'Variable  ',ibp,'  is fixed.'
      if (nleft .eq. 0 .and. nbreak .eq. n) then
c                                             all n variables are fixed,
c                                                return with xcp as GCP.
	 dtm = dt
	 goto 999
      endif
 
c     Update the derivative information.
 
      nint = nint + 1
      dibp2 = dibp**2
 
c     Update f1 and f2.
 
c        temporarily set f1 and f2 for col=0.
      f1 = f1 + dt*f2 + dibp2 - theta*dibp*zibp
      f2 = f2 - theta*dibp2

      if (col .gt. 0) then
c                          update c = c + dt*p.
	 call daxpy(col2,dt,p,1,c,1)
 
c           choose wbp,
c           the row of W corresponding to the breakpoint encountered.
      	 pointr = head
         do 70 j = 1,col
	    wbp(j) = wy(ibp,pointr)
	    wbp(col + j) = theta*ws(ibp,pointr)
            pointr = mod(pointr,m) + 1
  70     continue 
 
c           compute (wbp)Mc, (wbp)Mp, and (wbp)M(wbp)'.
         call bmv(m,sy,wt,col,wbp,v,info)
	 if (info .ne. 0) return
	 wmc = ddot(col2,c,1,v,1)
	 wmp = ddot(col2,p,1,v,1) 
	 wmw = ddot(col2,wbp,1,v,1)
 
c           update p = p - dibp*wbp. 
       	 call daxpy(col2,-dibp,wbp,1,p,1)
 
c           complete updating f1 and f2 while col > 0.
      	 f1 = f1 + dibp*wmc
      	 f2 = f2 + 2.0d0*dibp*wmp - dibp2*wmw
      endif

      if (nleft .gt. 0) then
         dtm = -f1/f2
         goto 777
c                 to repeat the loop for unsearched intervals. 
      else if(bnded) then
      	 f1 = zero
      	 f2 = zero
	 dtm = zero
      else
         dtm = -f1/f2
      endif 

c------------------- the end of the loop -------------------------------
 
 888  continue
      if (iprint .ge. 99) then
         write (6,*)
         write (6,*) 'GCP found in this segment'
         write (6,4010) nint,f1,f2
         write (6,6010) dtm
      endif 
      if (dtm .le. zero) dtm = zero
      tsum = tsum + dtm
 
c     Move free variables (i.e., the ones w/o breakpoints) and 
c       the variables whose breakpoints haven't been reached.
 
      call daxpy(n,tsum,d,1,xcp,1)
 
 999  continue
 
c     Update c = c + dtm*p = W'(x^c - x) 
c       which will be used in computing r = Z'(B(x^c - x) + g).
 
      if (col .gt. 0) call daxpy(col2,dtm,p,1,c,1)
      if (iprint .gt. 100) write (6,1010) (xcp(i),i = 1,n)
      if (iprint .ge. 99) write (6,2010)

 1010 format ('Cauchy X =  ',/,(4x,1p,6(1x,d11.4)))
 2010 format (/,'---------------- exit CAUCHY----------------------',/)
 3010 format (/,'---------------- CAUCHY entered-------------------')
 4010 format ('Piece    ',i3,' --f1, f2 at start point ',1p,2(1x,d11.4))
 4011 format (/,'Piece    ',i3,' --f1, f2 at start point ',
     +        1p,2(1x,d11.4))
 5010 format ('Distance to the next break point =  ',1p,d11.4)
 6010 format ('Distance to the stationary point =  ',1p,d11.4) 
 
      return
 
      end

c====================== The end of cauchy ==============================

      subroutine compr(n,m,x,f,g,ws,wy,sy,wt,z,r,d,wa,index,
     +                 theta,col,head,nfree,cnstnd,info)
 
      logical cnstnd
      integer n,m,col,head,nfree,info,index(n)
      double precision f,theta,
     +       x(n),g(n),z(n),r(n),d(n),wa(4*m),
     +       ws(n,m),wy(n,m),sy(m,m),wt(m,m)

c     ************
c
c     Subroutine compr 
c
c       This subroutine computes r=-Z'B(xcp-xk)-Z'g by using 
c         wa(2m+1)=W'(xcp-x) from subroutine cauchy.
c
c     Subprograms called:
c
c       L-BFGS-B Library ... bmv.
c
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************
 
      integer i,j,k,pointr
      double precision a1,a2

      if (.not. cnstnd .and. col .gt. 0) then 
         do 26 i = 1, n
	    r(i) = -g(i)
  26     continue
      else
         do 30 i = 1, nfree
            k = index(i)
	    r(i) = -theta*(z(k) - x(k)) - g(k)
  30     continue
     	 call bmv(m,sy,wt,col,wa(2*m+1),wa(1),info)
         if (info .ne. 0) then
            info = -8
	    return
         endif
     	 pointr = head 
     	 do 34 j = 1, col
       	    a1 = wa(j)
            a2 = theta*wa(col + j)
	    do 32 i = 1, nfree
	       k = index(i)
	       r(i) = r(i) + wy(k,pointr)*a1 + ws(k,pointr)*a2
  32        continue
	    pointr = mod(pointr,m) + 1
  34     continue
      endif

      return

      end

c======================= The end of compr ==============================

      subroutine errchk(n,m,factr,isbmin,l,u,nbd,task,info,k)
 
      character*60 task
      integer n,m,isbmin,info,k,nbd(n)
      double precision factr,l(n),u(n)

c     ************
c
c     Subroutine errchk
c
c     This subroutine checks the validity of the input data.
c
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************

      double precision one,zero
      parameter(one=1.0d0,zero=0.0d0)

      integer i

c     Check the input arguments for errors.

      if (n .le. 0) task = 'ERROR: N .LE. 0'
      if (m .le. 0) task = 'ERROR: M .LE. 0'
      if (factr .lt. zero) task = 'ERROR: FACTR .LT. 0'
      if (isbmin.lt.1 .or. isbmin.gt.3) task = 'ERROR: INVALID ISBMIN'

c     Check the validity of the arrays nbd(i), u(i), and l(i).

      do 10 i = 1, n
         if (nbd(i) .lt. 0 .or. nbd(i) .gt. 3) then
c                                                   return
            task = 'ERROR: INVALID NBD'
	    info = -6
	    k = i
         endif
	 if (nbd(i) .eq. 2) then
	    if (l(i) .gt. u(i)) then
c                                    return
               task = 'ERROR: NO FEASIBLE SOLUTION'
	       info = -7
	       k = i
            endif
         endif
  10  continue

      return

      end

c======================= The end of errchk =============================

      subroutine formh(n,m,sy,yy,sg,yg,sgo,ygo,rr,col,iupdat)
 
      integer n,m,col,iupdat
      double precision rr,sg(m),yg(m),sgo(m),ygo(m),sy(m,m),yy(m,m)


c     ************
c
c     Subroutine formh
c
c       This subroutine forms the middle matrix in H, update YY and
c         the upper triangular of SY.
c
c     Subprograms called:
c
c       Linpack ... dcopy.
c
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************

      integer i,j

c     Form the middle matrix in H, update YY,
c                                     and the upper triangular of SY.
 
      if (iupdat .gt. m) then
c                              move the old information.
         do 57 j = 1, col - 1
            call dcopy(col-1,yy(2,j+1),1,yy(1,j),1)
            call dcopy(j-1,sy(2,j+1),1,sy(1,j),1)
  57     continue
      endif
 
c     Add new information: last columns of YY and SY.
 
      do 58 i = 1, col - 1
         yy(i,col) = yg(i) - ygo(i)
	 sy(i,col) = sg(i) - sgo(i)
         yy(col,i) = yy(i,col)
  58  continue
      yy(col,col) = rr

      return

      end

c======================= The end of formh ==============================
 
      subroutine formk(n,nsub,ind,nenter,ileave,indx2,iupdat,
     +                 updatd,wn,wn1,m,ws,wy,sy,theta,col,head,info)

      integer n,nsub,m,col,head,nenter,ileave,iupdat,info,
     +        ind(n),indx2(n)
      double precision theta,
     +       wn(2*m,2*m),wn1(2*m,2*m),ws(n,m),wy(n,m),sy(m,m)
      logical updatd

c     ************
c
c     Subroutine formk 
c
c     This subroutine forms  the LEL^T factorization of the indefinite
c       matrix    K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
c                     [L_a -R_z           theta*S'AA'S ]
c                                                    where E = [-I  0]
c                                                              [ 0  I]
c     The matrix K can be shown to be equal to the matrix M^[-1]N
c       occurring in section 5.1 of the paper 
c       "A limited memory algorithm for bound constrained optimization",  
c        and to the matrix Mbar^[-1] Nbar in section 5.3.
c
c     The subroutine statement is:
c
c       subroutine formk(n,nsub,ind,nenter,ileave,indx2,iupdat,
c                        updatd,wn,wn1,m,ws,wy,sy,theta,col,head,info)
c
c       where
c
c     n is an integer variable.
c       On entry n is the dimension of the problem.
c       On exit n is unchanged.
c
c     nsub is an integer variable
c       On entry nsub is the number of subspace variables in free set.
c       On exit nsub is not changed.
c
c     ind is an integer array of dimension nsub.
c       On entry ind specifies the indices of subspace variables.
c       On exit ind is unchanged. 
c
c     nenter is an integer variable.
c       On entry nenter is the number of variables entering the 
c         free set.
c       On exit nenter is unchanged. 
c
c     ileave is an integer variable.
c       On entry indx2(ileave),...,indx2(n) are the variables leaving
c         the free set.
c       On exit ileave is unchanged. 
c
c     indx2 is an integer array of dimension n.
c       On entry indx2(1),...,indx2(nenter) are the variables entering
c         the free set, while indx2(ileave),...,indx2(n) are the
c         variables leaving the free set.
c       On exit indx2 is unchanged. 
c
c     iupdat is an integer variable.
c       On entry iupdat is the total number of BFGS updates made so far.
c       On exit iupdat is unchanged. 
c
c     updatd is a logical variable.
c       On entry 'updatd' is true if the L-BFGS matrix is updatd.
c       On exit 'updatd' is unchanged. 
c
c     wn is a double precision array of dimension 2m x 2m.
c       On entry wn is unspecified.
c       On exit the upper triangle of wn stores the LEL^T factorization
c         of the 2*col x 2*col indefinite matrix
c                     [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
c                     [L_a -R_z           theta*S'AA'S ]
c
c     wn1 is a double precision array of dimension 2m x 2m.
c       On entry wn1 stores the lower triangular part of 
c                     [Y' ZZ'Y   L_a'+R_z']
c                     [L_a+R_z   S'AA'S   ]
c         in the previous iteration.
c       On exit wn1 stores the corresponding updated matrices.
c       The purpose of wn1 is just to store these inner products
c       so they can be easily updated and inserted into wn.
c
c     m is an integer variable.
c       On entry m is the maximum number of variable metric corrections
c         used to define the limited memory matrix.
c       On exit m is unchanged.
c
c     ws, wy, sy, and wtyy are double precision arrays;
c     theta is a double precision variable;
c     col is an integer variable;
c     head is an integer variable.
c       On entry they store the information defining the
c                                          limited memory BFGS matrix:
c         ws(n,m) stores S, a set of s-vectors;
c         wy(n,m) stores Y, a set of y-vectors;
c         sy(m,m) stores S'Y;
c         wtyy(m,m) stores the cholesky factorization
c                                   of (theta*S'S+LD^(-1)L')
c         theta is the scaling factor specifying B_0 = theta I;
c         col is the number of variable metric corrections stored;
c         head is the location of the 1st s- (or y-) vector in S (or Y).
c       On exit they are unchanged.
c
c     info is an integer variable.
c       On entry info is unspecified.
c       On exit info =  0 for normal return;
c                    = -1 when the 1st Cholesky factorization failed;
c                    = -2 when the 2st Cholesky factorization failed.
c
c     Subprograms called:
c
c       Linpack ... dcopy, dpofa, dtrsl.
c
c
c     References:
c       R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
c       memory algorithm for bound constrained optimization'',
c       SIAM J. Scientific Computing 16 (1995), no. 5.
c
c       C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a
c       limited memory FORTRAN code for solving bound constrained
c       optimization problems'', Tech. Report, EECS Department,
c       Northwestern University, 1994.
c
c       (Postscript files of these papers are available via anonymous
c        ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************

      double precision one,zero
      parameter(one=1.0d0,zero=0.0d0)

      integer m2,ipntr,jpntr,iy,is,jy,js,is1,js1,k1,i,k,col2,
     +        pbegin,pend,dbegin,dend,upcl
      double precision ddot,temp1,temp2,temp3,temp4
 
c     Form the lower triangular part of
c               WN1 = [Y' ZZ'Y   L_a'+R_z'] 
c                     [L_a+R_z   S'AA'S   ]
c        where L_a is the strictly lower triangular part of S'AA'Y
c              R_z is the upper triangular part of S'ZZ'Y.
      
      if (updatd) then
         if (iupdat .gt. m) then 
c                                 shift old part of WN1.
            do 10 jy = 1, m - 1
               js = m + jy
	       call dcopy(m-jy,wn1(jy+1,jy+1),1,wn1(jy,jy),1)
 	       call dcopy(m-jy,wn1(js+1,js+1),1,wn1(js,js),1)
 	       call dcopy(m-1,wn1(m+2,jy+1),1,wn1(m+1,jy),1)
  10        continue
         endif
 
c          put new rows in blocks (1,1), (2,1) and (2,2).
         pbegin = 1
	 pend = nsub
         dbegin = nsub + 1
	 dend = n
         iy = col
         is = m + col
         ipntr = head + col - 1
         if (ipntr .gt. m) ipntr = ipntr - m	
         jpntr = head
         do 20 jy = 1, col
            js = m + jy
            temp1 = zero
	    temp2 = zero
	    temp3 = zero
c             compute element jy of row 'col' of Y'ZZ'Y
	    do 15 k = pbegin, pend
	       k1 = ind(k)
	       temp1 = temp1 + wy(k1,ipntr)*wy(k1,jpntr)
  15        continue
c             compute elements jy of row 'col' of L_a and S'AA'S
	    do 16 k = dbegin, dend
	       k1 = ind(k)
	       temp2 = temp2 + ws(k1,ipntr)*ws(k1,jpntr)
	       temp3 = temp3 + ws(k1,ipntr)*wy(k1,jpntr)
  16        continue
	    wn1(iy,jy) = temp1
	    wn1(is,js) = temp2
	    wn1(is,jy) = temp3
            jpntr = mod(jpntr,m) + 1
  20     continue
 
c          put new column in block (2,1).
         jy = col	
         jpntr = head + col - 1
         if (jpntr .gt. m) jpntr = jpntr - m
         ipntr = head
         do 30 i = 1, col
            is = m + i
	    temp3 = zero
c             compute element i of column 'col' of R_z
	    do 25 k = pbegin, pend
	       k1 = ind(k)
	       temp3 = temp3 + ws(k1,ipntr)*wy(k1,jpntr)
  25        continue 
	    ipntr = mod(ipntr,m) + 1
            wn1(is,jy) = temp3
  30     continue
	 upcl = col - 1
      else
         upcl = col
      endif
 
c       modify the old parts in blocks (1,1) and (2,2) due to changes
c       in the set of free variables.
      ipntr = head	
      do 45 iy = 1, upcl
         is = m + iy
	 jpntr = head
      	 do 40 jy = 1, iy
	    js = m + jy
	    temp1 = zero
	    temp2 = zero
	    temp3 = zero
	    temp4 = zero
	    do 35 k = 1, nenter
	       k1 = indx2(k)
	       temp1 = temp1 + wy(k1,ipntr)*wy(k1,jpntr)
	       temp2 = temp2 + ws(k1,ipntr)*ws(k1,jpntr)
  35        continue
	    do 36 k = ileave, n
	       k1 = indx2(k)
	       temp3 = temp3 + wy(k1,ipntr)*wy(k1,jpntr)
	       temp4 = temp4 + ws(k1,ipntr)*ws(k1,jpntr)
  36        continue
	    wn1(iy,jy) = wn1(iy,jy) + temp1 - temp3 
	    wn1(is,js) = wn1(is,js) - temp2 + temp4 
	    jpntr = mod(jpntr,m) + 1
  40     continue
         ipntr = mod(ipntr,m) + 1
  45  continue
 
c       modify the old parts in block (2,1).
      ipntr = head      
      do 60 is = m + 1, m + upcl
         jpntr = head 
         do 55 jy = 1, upcl
            temp1 = zero
	    temp3 = zero
	    do 50 k = 1, nenter
	       k1 = indx2(k)
	       temp1 = temp1 + ws(k1,ipntr)*wy(k1,jpntr)
  50	    continue
	    do 51 k = ileave, n
	       k1 = indx2(k)
	       temp3 = temp3 + ws(k1,ipntr)*wy(k1,jpntr)
  51	    continue
         if (is .le. jy + m) then
	       wn1(is,jy) = wn1(is,jy) + temp1 - temp3  
	    else
	       wn1(is,jy) = wn1(is,jy) - temp1 + temp3  
	    endif
	    jpntr = mod(jpntr,m) + 1
  55     continue
         ipntr = mod(ipntr,m) + 1
  60  continue
 
c     Form the upper triangle of WN = [D+Y' ZZ'Y/theta   -L_a'+R_z' ] 
c                                     [-L_a +R_z        S'AA'S*theta]

      m2 = 2*m
      do 70 iy = 1, col
	 is = col + iy
	 is1 = m + iy
      	 do 65 jy = 1, iy
	    js = col + jy
            js1 = m + jy
 	    wn(jy,iy) = wn1(iy,jy)/theta
 	    wn(js,is) = wn1(is1,js1)*theta
  65     continue
      	 do 66 jy = 1, iy - 1
 	    wn(jy,is) = -wn1(is1,jy)
  66     continue
      	 do 67 jy = iy, col
 	    wn(jy,is) = wn1(is1,jy)
  67     continue
 	 wn(iy,iy) = wn(iy,iy) + sy(iy,iy)
  70  continue

c     Form the upper triangle of WN= [  LL'            L^-1(-L_a'+R_z')] 
c                                    [(-L_a +R_z)L'^-1   S'AA'S*theta  ]

c        first Cholesky factorization of (1,1) block of wn to get LL'
c                          with L' stored in the upper triangle of wn.
      call dpofa(wn,m2,col,info)
      if (info .ne. 0) then
	 info = -1
	 return
      endif
c        then form L^-1(-L_a'+R_z') in the (1,2) block.
      col2 = 2*col
      do 71 js = col+1 ,col2
         call dtrsl(wn,m2,col,wn(1,js),11,info)
  71  continue

c     Form S'AA'S*theta + (L^-1(-L_a'+R_z'))'L^-1(-L_a'+R_z') in the
c        upper triangle of (2,2) block of wn.
                      

      do 72 is = col+1, col2
         do 74 js = is, col2
	       wn(is,js) = wn(is,js) + ddot(col,wn(1,is),1,wn(1,js),1)
  74        continue
  72     continue

c     Cholesky factorization of (2,2) block of wn.

      call dpofa(wn(col+1,col+1),m2,col,info)
      if (info .ne. 0) then
	 info = -2
	 return
      endif

      return

      end

c======================= The end of formk ==============================

      subroutine formt(n,m,wt,sy,ss,col,theta,info)
 
      integer n,m,col,info
      double precision theta,wt(m,m),sy(m,m),ss(m,m)

c     ************
c
c     Subroutine formt
c
c       This subroutine forms the upper half of the pds
c         T = theta*SS + L*D^(-1)*L', stores T in the upper triangular
c         of the array wt, and performs the Cholesky factorization of T
c         to produce J*J', with J' stored in the upper triangular of wt.
c
c     Subprograms called:
c
c       Linpack ... dpofa.
c
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************

      double precision zero
      parameter(zero=0.0d0)

      integer i,j,k,k1
      double precision ddum

c     Form the upper half of the pds T = theta*SS + L*D^(-1)*L',
c        store T in the upper triangular of the array wt.
 
      do 52 j = 1, col
      	 wt(1,j) = theta*ss(1,j)
  52  continue
      do 55 i = 2, col
	 do 54 j = i, col
            k1 = min(i,j) - 1
            ddum  = zero
            do 53 k = 1, k1
               ddum  = ddum + sy(i,k)*sy(j,k)/sy(k,k)
  53        continue
            wt(i,j) = ddum + theta*ss(i,j)
  54     continue
  55  continue
 
c     Cholesky factorize T to J*J' with 
c        J' stored in the upper triangular of wt.
 
      call dpofa(wt,m,col,info)
      if (info .ne. 0) then
         info = -3
      endif

      return

      end

c======================= The end of formt ==============================
 
      subroutine freev(n,nfree,index,nenter,ileave,indx2,
     +                 iwhere,wrk,updatd,cnstnd,iprint,iter)

      integer n,nfree,nenter,ileave,iprint,iter,
     +        index(n),indx2(n),iwhere(n)
      logical wrk,updatd,cnstnd

c     ************
c
c     Subroutine freev 
c
c     This subroutine counts the entering and leaving variables when
c       iter > 0, and finds the index set of free and active variables
c       at the GCP.
c
c     cnstnd is a logical variable indicating whether bounds are present
c
c     index is an integer array of dimension n
c       for i=1,...,nfree, index(i) are the indices of free variables
c       for i=nfree+1,...,n, index(i) are the indices of bound variables
c       On entry after the first iteration, index gives 
c         the free variables at the previous iteration.
c       On exit it gives the free variables based on the determination
c         in cauchy using the array iwhere.
c
c     indx2 is an integer array of dimension n
c       On entry indx2 is unspecified.
c       On exit with iter>0, indx2 indicates which variables
c          have changed status since the previous iteration.
c       For i= 1,...,nenter, indx2(i) have changed from bound to free.
c       For i= ileave+1,...,n, indx2(i) have changed from free to bound.
c 
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************
 
      integer iact,i,k

      nenter = 0
      ileave = n + 1
      if (iter .gt. 0 .and. cnstnd) then
c                           count the entering and leaving variables.
	 do 20 i = 1, nfree
	    k = index(i)
	    if (iwhere(k) .gt. 0) then
	       ileave = ileave - 1
	       indx2(ileave) = k
	       if (iprint .ge. 100) write (6,*)
     +             'Variable ',k,' leaves the set of free variables'
            endif
  20     continue
	 do 22 i = 1 + nfree, n
	    k = index(i)
	    if (iwhere(k) .le. 0) then
	       nenter = nenter + 1
	       indx2(nenter) = k
	       if (iprint .ge. 100) write (6,*)
     +             'Variable ',k,' enters the set of free variables'
            endif
  22     continue
         if (iprint .ge. 99) write (6,*)
     +       n+1-ileave,' variables leave; ',nenter,' variables enter'
      endif
      wrk = (ileave .lt. n+1) .or. (nenter .gt. 0) .or. updatd
 
c     Find the index set of free and active variables at the GCP.
 
      nfree = 0 
      iact = n + 1
      do 24 i = 1, n
	 if (iwhere(i) .le. 0) then
	    nfree = nfree + 1
	    index(nfree) = i
         else
            iact = iact - 1
            index(iact) = i
         endif
  24  continue
      if (iprint .ge. 99) write (6,*)
     +      nfree,' variables are free at GCP ',iter + 1  

      return

      end

c======================= The end of freev ==============================

      subroutine hmv(m,sy,yy,theta,col,v,p,info)
 
      integer m,col,info
      double precision theta,
     +                 sy(m,m),yy(m,m),v(2*col),p(2*col)

c     ************
c
c     Subroutine hmv
c
c     This subroutine computes the product of the 2m x 2m center matrix.
c       in the compact L-BFGS formula of H with a 2m vector v.  
c       It returns the product in p.
c
c     The subroutine statement is:
c
c       subroutine hmv(m,sy,yy,theta,col,v,p)
c
c       where
c
c     m is an integer variable.
c       On entry m is the maximum number of variable metric corrections
c         used to define the limited memory matrix.
c       On exit m is unchanged.
c
c     sy is a double precision array of dimension m x m.
c       On entry sy specifies the matrix S'Y.
c       On exit sy is unchanged.
c
c     yy is a double precision array of dimension m x m.
c       On entry yy specifies the matrix Y'Y.
c       On exit yy is unchanged.
c
c     col is an integer variabl.
c       On entry col specifies the number of s-vectors (or y-vectors)
c         stored in the compact L-BFGS formula.
c       On exit col is unchanged.
c
c     theta is a double precision variable.
c       On entry theta is the scaling factor specifying H_0 = I/theta.
c       On exit theta is unchanged.
c
c     v is a double precision array of dimension 2col.
c       On entry v specifies vector v.
c       On exit v is unchanged.
c
c     p is a double precision array of dimension 2col.
c       On entry p is unspecified.
c       On exit p is the product Mv.
c
c     info is an integer variable.
c       On entry info is unspecified.
c       On exit info = 0       for normal return,
c                    = nonzero for abnormal return when the system
c                                     in call of dtrsl is singular.
c
c     Subprograms called:
c
c       Linpack ... dtrsl.
c
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************

      integer i,i2
      double precision temp,ddot
c
      if (col .eq. 0) return
c
c     Compute [ p1 ] = [ 0         -R^(-1)                 ] [ v1 ]
c             [ p2 ]   [ -R^(-T)   R^(-T)(D+YY/theta)R^(-1)] [ v2 ] 
 
c       p1:=R^(-1)(-v2); p2:=-v1
      do 10 i = 1, col
         i2 = col + i
         temp = -v(i)
         p(i) = -v(i2)
  	 p(i2) = temp
  10  continue
      call dtrsl(sy,m,col,p,01,info)
      if (info .ne. 0) return
 
c       p2:=R'^(-1)(p2-(D+YY/theta)p1)
      do 20 i = 1, col
         i2 = col + i
         temp = ddot(col,yy(i,1),m,p,1)
         p(i2) = p(i2) - sy(i,i)*p(i) - temp/theta
  20   continue  
      call dtrsl(sy,m,col,p(col+1),11,info)
      if (info .ne. 0) return

      return

      end

c======================== The end of hmv ===============================

      subroutine lnsrch(n,l,u,nbd,x,f,fold,gd,gdold,g,d,r,t,z,stp,
     +                  dnorm,dtd,xstep,stpmx,iter,ifun,iback,nfgv,
     +                  info,task,boxed,cnstnd,csave,isave,dsave)

      character*60 task,csave
      logical boxed,cnstnd
      integer n,iter,ifun,iback,nfgv,info,nbd(n),isave(2)
      double precision f,fold,gd,gdold,stp,dnorm,dtd,xstep,stpmx,
     +                 x(n),l(n),u(n),g(n),d(n),r(n),t(n),z(n),dsave(13)
c     **********
c
c     Subroutine lnsrch
c
c     This subroutine performs line search within the feasible region 
c       in the direction specified by d.
c
c     Subprograms called:
c
c       Minpack2 Library ... dcsrch.
c
c       Linpack ... dtrsl, ddot.
c
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     **********

      double precision one,zero,big
      parameter(one=1.0d0,zero=0.0d0,big=1.0d+10)
      double precision ftol,gtol,xtol
      parameter(ftol=1.0d-3,gtol=0.9d0,xtol=0.1d0)

      integer i
      double precision ddot,a1,a2

      if (task(1:5) .eq. 'FG_LN') goto 556

      dtd = ddot(n,d,1,d,1)
      dnorm = sqrt(dtd)

c     Determine the maximum step length.

      stpmx = big
      if (cnstnd) then
         if (iter .eq. 0) then
            stpmx = one
         else
            do 43 i = 1, n
               a1 = d(i)
               if (nbd(i) .ne. 0) then
                  if (a1 .lt. zero .and. nbd(i) .le. 2) then
                     a2 = l(i) - x(i)
                     if (a2 .ge. zero) then
                        stpmx = zero
                     else if (a1*stpmx .lt. a2) then
                        stpmx = a2/a1
                     endif
                  else if (a1 .gt. zero .and. nbd(i) .ge. 2) then
                     a2 = u(i) - x(i)
                     if (a2 .le. zero) then
                        stpmx = zero
                     else if (a1*stpmx .gt. a2) then
                        stpmx = a2/a1
                     endif
                  endif
               endif
  43        continue
         endif
      endif
 
      if (iter .eq. 0 .and. .not. boxed) then
         stp = min(one/dnorm, stpmx)
      else
	 stp = one
      endif 

      call dcopy(n,x,1,t,1)
      call dcopy(n,g,1,r,1)
      fold = f
      ifun = 0
      iback = 0
      csave = 'START'
 556  continue
      gd = ddot(n,g,1,d,1)
      if (ifun .eq. 0) then
	 gdold=gd
         if (gd .ge. zero) then
c                               the directional derivative >=0.
c                               Line search is impossible.
            info = -4
            return
         endif
      endif

      call dcsrch(f,gd,stp,ftol,gtol,xtol,zero,stpmx,csave,isave,dsave)

      xstep = stp*dnorm
      if (csave(1:4) .ne. 'CONV' .and. csave(1:4) .ne. 'WARN') then
	 task = 'FG_LNSRCH'
	 ifun = ifun + 1
         nfgv = nfgv + 1
         iback = ifun - 1 
         if (stp .eq. one) then
            call dcopy(n,z,1,x,1)
         else
            do 41 i = 1, n
	       x(i) = stp*d(i) + t(i)
  41        continue
         endif
c          return to the driver for calculating f and g.
	 return
      else
         task = 'NEW_X'
      endif

      return

      end

c======================= The end of lnsrch =============================

      subroutine matupd(n,m,ws,wy,sy,ss,d,r,itail,
     +                  iupdat,col,head,theta,rr,dr,stp,dtd)
 
      integer n,m,itail,iupdat,col,head
      double precision theta,rr,dr,stp,dtd,d(n),r(n),
     +       ws(n,m),wy(n,m),sy(m,m),ss(m,m)

c     ************
c
c     Subroutine matupd
c
c       This subroutine updates matrices WS and WY, and form the
c         middle matrix in B.
c
c     Subprograms called:
c
c       Linpack ... dcopy, ddot.
c
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************
 
      integer j,pointr
      double precision ddot
      double precision one
      parameter(one=1.0d0)

c     Set pointers for matrices WS and WY.
 
      if (iupdat .le. m) then
	 col = iupdat
	 itail = mod(head+iupdat-2,m) + 1
      else
	 itail = mod(itail,m) + 1
	 head = mod(head,m) + 1
      endif
 
c     Update matrices WS and WY.

      call dcopy(n,d,1,ws(1,itail),1)
      call dcopy(n,r,1,wy(1,itail),1)
 
c     Set theta=yy/ys.
 
      theta = rr/dr
 
c     Form the middle matrix in B.
 
c        update the upper triangular of SS,
c                                         and the lower triagular of SY:
      if (iupdat .gt. m) then
c                              move old information
         do 50 j = 1, col - 1
            call dcopy(j,ss(2,j+1),1,ss(1,j),1)
            call dcopy(col-j,sy(j+1,j+1),1,sy(j,j),1)
  50     continue
      endif
c        add new information: the last row of SY
c                                             and the last column of SS:
      pointr = head
      do 51 j = 1, col - 1
	 sy(col,j) = ddot(n,d,1,wy(1,pointr),1)
	 ss(j,col) = ddot(n,ws(1,pointr),1,d,1)
         pointr = mod(pointr,m) + 1
  51  continue
      if (stp .eq. one) then
         ss(col,col) = dtd
      else
         ss(col,col) = stp*stp*dtd
      endif
      sy(col,col) = dr
 
      return

      end

c======================= The end of matupd =============================

      subroutine print1(n,m,l,u,x,iprint,isbmin,itfile,epsmch)
 
      integer n,m,iprint,isbmin,itfile
      double precision epsmch,x(n),l(n),u(n)

c     ************
c
c     Subroutine print1
c
c     This subroutine print the input data, initial point, upper and
c       lower bounds of each variable, machine precision, as well as 
c       the headings of the output.
c
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************

      integer i

      if (iprint .ge. 0) then
         write (6,7001) epsmch
         if (isbmin .eq. 2) then
	    write (6,8002)
	 else if (isbmin .eq. 1) then
	    write (6,8001)
	 else
	    write (6,8003)
	 endif
         write (6,*) 'N = ',n,'    M = ',m
         if (iprint .ge. 1) then
            write (itfile,2001) epsmch
            if (isbmin .eq. 2) then
	       write (itfile,8002)
	    else if (isbmin .eq. 1) then
	       write (itfile,8001)
	    else
	       write (itfile,8003)
	    endif
            write (itfile,*)'N = ',n,'    M = ',m
	    write (itfile,9001)
            if (iprint .gt. 100) then
               write (6,1004) 'L =',(l(i),i = 1,n)
               write (6,1004) 'X0 =',(x(i),i = 1,n)
               write (6,1004) 'U =',(u(i),i = 1,n)
            endif 
         endif
      endif 

 1004 format (/,a4, 1p, 6(1x,d11.4),/,(4x,1p,6(1x,d11.4)))
 2001 format ('RUNNING THE L-BFGS-B CODE',/,/,
     + 'it    = iteration number',/,
     + 'nf    = number of function evaluations',/,
     + 'nint  = number of segments explored during the Cauchy search',/,
     + 'nact  = number of active bounds at the generalized Cauchy point'
     + ,/,
     + 'sub   = manner in which the subspace minimization terminated:'
     + ,/,'        con = converged, bnd = a bound was reached',/,
     + 'itls  = number of iterations performed in the line search',/,
     + 'stepl = step length used',/,
     + 'tstep = norm of the displacement (total step)',/,
     + 'projg = norm of the projected gradient',/,
     + 'f     = function value',/,/,
     + '           * * *',/,/,
     + 'Machine precision =',1p,d10.3)
 7001 format ('RUNNING THE L-BFGS-B CODE',/,/,
     + '           * * *',/,/,
     + 'Machine precision =',1p,d10.3)
 8001 format ('Subspace minimization done by Direct method',/)
 8002 format ('Subspace minimization done by Dual method',/)
 8003 format ('Subspace minimization done by CG',/)
 9001 format (/,3x,'it',3x,'nf',2x,'nint',2x,'nact',2x,'sub',2x,'itls',
     +        2x,'stepl',4x,'tstep',5x,'projg',8x,'f')

      return

      end

c======================= The end of print1 =============================

      subroutine print2(n,x,f,g,iprint,itfile,iter,nfgv,nact,
     +                  sbgnrm,nint,word,iword,iback,stp,xstep)
 
      character*3 word
      integer n,iprint,itfile,iter,nfgv,nact,nint,iword,iback
      double precision f,sbgnrm,stp,xstep,x(n),g(n)

c     ************
c
c     Subroutine print2
c
c     This subroutine print out new information after a succesful
c       line search. 
c
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************

      integer i,imod

c           'word' records the status of subspace solutions.
      if (iword .eq. 0) then
c                            the subspace minimizatin converges.
	 word = 'con'
      else if (iword .eq. 1) then
c                             the subspace minimisatin stops at a bound.
         word = 'bnd'
      else if (iword .eq. 5) then
c                                 the truncated Newton has been used.
	 word = 'TNT'
      else
         word = '---'
      endif
      if (iprint .ge. 99) then
         write (6,*) 'LINE SEARCH',iback,' times; norm of step = ',xstep
         write (6,2001) iter,f,sbgnrm
         if (iprint .gt. 100) then	
            write (6,1004) 'X =',(x(i), i = 1, n)
            write (6,1004) 'G =',(g(i), i = 1, n)
         endif
      else if (iprint .gt. 0) then 
         imod = mod(iter,iprint)
         if (imod .eq. 0) write (6,2001) iter,f,sbgnrm
      endif
      if (iprint .ge. 1) write (itfile,3001)
     +          iter,nfgv,nint,nact,word,iback,stp,xstep,sbgnrm,f

 1004 format (/,a4, 1p, 6(1x,d11.4),/,(4x,1p,6(1x,d11.4)))
 2001 format
     +  (/,'At iterate',i5,4x,'f= ',1p,d12.5,4x,'|proj g|= ',1p,d12.5)
 3001 format(2(1x,i4),2(1x,i5),2x,a3,1x,i4,1p,2(2x,d7.1),1p,2(1x,d10.3))

      return

      end

c======================= The end of print2 =============================

      subroutine print3(n,l,u,nbd,x,f,task,iprint,info,itfile,
     +                  iter,nfgv,nintol,ncgtol,nskip,nact,sbgnrm,
     +                  time,nint,word,iback,stp,xstep,k,
     +                  cachyt,sbtime,lnscht)
 
      character*60 task
      character*3 word
      integer n,iprint,info,itfile,iter,nfgv,nintol,ncgtol,nskip,
     +        nact,nint,iback,k,nbd(n)
      double precision f,sbgnrm,time,stp,xstep,cachyt,sbtime,lnscht,
     +                 x(n),l(n),u(n)

c     ************
c
c     Subroutine print3
c
c     This subroutine print out the information when either some inbuilt
c       convergence criteria are satisfied or some error message is
c       generated.
c       
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************

      integer i

      if (task(1:5) .eq. 'ERROR') goto 999

      if (iprint .ge. 0) then
         write (6,3003)
         write (6,3004)
         write(6,3005) n,iter,nfgv,nintol,ncgtol,nskip,nact,sbgnrm,f
         if (iprint .ge. 100) then
            write (6,1004) 'X =',(x(i),i = 1,n)
         endif  
         if (iprint .ge. 1) write (6,*) ' F =',f
      endif 
 999  continue
      if (iprint .ge. 0) then
         write (6,3009) task
         if (info .ne. 0) then
            if (info .eq. -1) write (6,9011)
            if (info .eq. -2) write (6,9012)
            if (info .eq. -3) write (6,9013)
            if (info .eq. -4) write (6,9014)
            if (info .eq. -5) write (6,9015)
            if (info .eq. -6) write (6,*)' Input nbd(',k,') is invalid.'
            if (info .eq. -7) 
     +	       write (6,*)' l(',k,') > u(',k,').  No feasible solution.'
            if (info .eq. -8) write (6,9018)
            if (info .eq. -9) write (6,9019)
         endif
         if (iprint .ge. 1) write (6,3007) cachyt,sbtime,lnscht
         write (6,3008) time
         if (iprint .ge. 1) then
            if (info .eq. -4 .or. info .eq. -9) then
               write (itfile,3002)
     +             iter,nfgv,nint,nact,word,iback,stp,xstep
            endif
            write (itfile,3009) task
            if (info .ne. 0) then
               if (info .eq. -1) write (itfile,9011)
               if (info .eq. -2) write (itfile,9012)
               if (info .eq. -3) write (itfile,9013)
               if (info .eq. -4) write (itfile,9014)
               if (info .eq. -5) write (itfile,9015)
               if (info .eq. -8) write (itfile,9018)
               if (info .eq. -9) write (itfile,9019)
            endif
            write (itfile,3008) time
         endif
      endif

 1004 format (/,a4, 1p, 6(1x,d11.4),/,(4x,1p,6(1x,d11.4)))
 3002 format(2(1x,i4),2(1x,i5),2x,a3,1x,i4,1p,2(2x,d7.1),6x,'-',10x,'-')
 3003 format (/,
     + '           * * *',/,/,
     + 'Tit   = total number of iterations',/,
     + 'Tnf   = total number of function evaluations',/,
     + 'Tnint = total number of segments explored during',
     +           ' Cauchy searches',/,
     + 'Titcg = Total number of CG iterations',/,
     + 'Skip  = number of BFGS updates skipped',/,
     + 'Nact  = number of active bounds at final generalized',
     +          ' Cauchy point',/,
     + 'Projg = norm of the final projected gradient',/,
     + 'F     = final function value',/,/,
     + '           * * *')
 3004 format (/,3x,'N',3x,'Tit',2x,'Tnf',2x,'Tnint',2x,'Titcg',1x,
     +       'Skip',2x,'Nact',5x,'Projg',8x,'F')
 3005 format (i5,2(1x,i4),2(1x,i6),(1x,i4),(1x,i5),1p,2(2x,d10.3))
 3006 format (i5,2(1x,i4),2(1x,i6),(1x,i4),(1x,i5),7x,'-',10x,'-')
 3007 format (/,' Cauchy                time',1p,e10.3,' seconds.',/ 
     +        ' Subspace minimization time',1p,e10.3,' seconds.',/
     +        ' Line search           time',1p,e10.3,' seconds.')
 3008 format (/,' Total User time',1p,e10.3,' seconds.',/)
 3009 format (/,a60)
 9011 format (/,
     +' Matrix in 1st Cholesky factorization in formk is not Pos. Def.')
 9012 format (/,
     +' Matrix in 2st Cholesky factorization in formk is not Pos. Def.')
 9013 format (/,
     +' Matrix in the Cholesky factorization in formt is not Pos. Def.')
 9014 format (/,
     +' Derivative >= 0, backtracking line search impossible.',/,
     +'   Previous x, f and g restored.',/,
     +' Possible causes: 1 error in function or gradient evaluation;',/,
     +'                  2 rounding error blurs the real information.')
 9015 format (/,
     +' Warning: there are more than 10 function and gradient',/,
     +'   evaluations in the last line search.  The termination',/,
     +'   may possibly be caused by a bad search direction.')
 9018 format (/,' The triangular system is singular.')
 9019 format (/,
     +' Line search cannot locate an adequate point after 20 function',/
     +,'  and gradient evaluations.  Previous x, f and g restored.',/,
     +' Possible causes: 1 error in function or gradient evaluation;',/,
     +'                  2 rounding error blurs the real information.')

      return

      end

c======================= The end of print3 =============================

      subroutine projgr(n,l,u,nbd,x,g,sbgnrm)

      integer n,nbd(n)
      double precision sbgnrm,x(n),l(n),u(n),g(n)

c     ************
c
c     Subroutine projgr
c
c     This subroutine computes the infinity norm of the projected
c       gradient.
c
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************

      double precision one,zero
      parameter(one=1.0d0,zero=0.0d0)

      integer i
      double precision gi

      sbgnrm = zero
      do 15 i = 1, n
	gi = g(i)
        if (nbd(i) .ne. 0) then
           if (gi .lt. zero) then
              if (nbd(i) .ge. 2) gi = max((x(i)-u(i)),gi)
       	   else
              if (nbd(i) .le. 2) gi = min((x(i)-l(i)),gi)
           endif
        endif
	sbgnrm = max(sbgnrm,abs(gi))
  15  continue

      return

      end

c======================= The end of projgr =============================
 
      subroutine subbv(n,m,nsub,ind,ws,wy,sy,wt,
     +                                    theta,col,head,v,p,wa,info) 

      integer n,m,nsub,col,head,info,
     +        ind(nsub)
      double precision theta,
     +       ws(n,m),wy(n,m),sy(m,m),wt(m,m),v(nsub),p(nsub),wa(2*col)

c     ************
c
c     Subroutine subbv 
c
c     Given a subspace specified by its size 'nsub' and an index set
c       'ind', as well as a vector v on this subspace, and the L-BFGS
c       Hessian approximation B, in terms of  ws, wy, sy, wt, theta,
c       col and head, this subroutine computes the product p = Z'BZv.
c 	
c     The subroutine statement is:
c
c       subroutine subbv(n,m,nsub,ind,ws,wy,sy,wt,
c    +                                    theta,col,head,v,p,wa,info) 
c
c       where
c
c     n is an integer variable.
c       On entry n is the dimension of the problem.
c       On exit n is unchanged.
c
c     m is an integer variable.
c       On entry m is the maximum number of variable metric corrections
c         used to define the limited memory matrix.
c       On exit m is unchanged.
c
c     nsub is an integer variable
c       On entry nsub specifies the size of the subspace.
c       On exit nsub is unchanged.
c
c     ind is an integer array of dimension nsub.
c       On entry ind specifies the indices of subspace variables.
c       On exit ind is unchanged.
c
c     v is a double precision array of dimension nsub.
c       On entry v specifies vector v.
c       On exit v is unchanged.
c
c     ws, wy, sy, and wt are double precision arrays;
c     theta is a double precision variable;
c     col is an integer variable;
c     head is an integer variable.
c       On entry they store the information defining the
c                                          limited memory BFGS matrix:
c         ws(n,m) stores S, a set of s-vectors;
c         wy(n,m) stores Y, a set of y-vectors;
c         sy(m,m) stores S'Y;
c         wt(m,m) stores the chelosky factorization
c                                             of (theta*S'S+LD^(-1)L').
c         theta is the scaling factor specifying B_0 = theta I;
c         col is the number of variable metric corrections stored;
c         head is the location of the 1st s- (or y-) vector in S (or Y).
c       On exit they are unchanged.
c
c     v is a double precision working array of dimension nsub.
c       On entry v is the vector v.
c       On exit v is unchanged.
c
c     p is a double precision working array of dimension nsub.
c       On entry p is unspecified.
c       On exit p is the product Z^(T)BZ*v.
c
c     wa is a double precision working array of dimension 2col.
c
c     info is an integer variable.
c       On entry info is 0.
c       On exit info = 0       for normal return,
c                    = nonzero for abnormal return when the system
c                                      in call of bmv is singular.
c
c     Subprograms called:
c
c       L-BFGS-B Library ... bmv.
c
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************
 
      double precision one,zero
      parameter(one=1.0d0,zero=0.0d0)
      
      integer i,j,k,pointr
      double precision sum1,sum2,a1,a2
 
      if (col .eq. 0) return
 
c     Compute wa = [      Y'Z*v ]
c                  [theta*S'Z*v ]
 
      pointr = head 
      do 20 i = 1, col
      	 sum1 = zero
	 sum2 = zero
 	 do 10 j = 1, nsub
 	    k = ind(j)
 	    sum1 = sum1 + wy(k,pointr)*v(j)
            sum2 = sum2 + ws(k,pointr)*v(j)
  10     continue
	 wa(i) = sum1
	 wa(col+i) = sum2*theta
         pointr = mod(pointr,m) + 1
  20  continue
 
c     Compute wa = M*wa.
 
      call bmv(m,sy,wt,col,wa,wa,info) 
      if (info .ne. 0) return
 
c     Compute p = theta*v-[Z'Y  theta*Z'S]*wa.
 
      do 30 i = 1, nsub
	 p(i) = theta*v(i)
  30  continue
      pointr = head 
      do 50 j = 1, col
         a1 = wa(j)
         a2 = theta*wa(col+j)
	 do 40 i = 1, nsub
	    k = ind(i)
	    p(i) = p(i) - wy(k,pointr)*a1 - ws(k,pointr)*a2
  40  	 continue
         pointr = mod(pointr,m) + 1
  50  continue

      return

      end

c====================== The end of subbv ===============================

      subroutine subcg(n,nsub,ind,l,u,nbd,x,r,p,q,m,ws,wy,sy,wt,
     +                 theta,col,head,itercg,wv,iword,iprint,info)
                                                                      
      integer n,nsub,m,col,head,itercg,iword,iprint,info,
     +        ind(nsub),nbd(n)
      double precision theta,
     +       l(n),u(n),x(n),r(n),p(nsub),q(nsub),wv(2*m),
     +       ws(n,col),wy(n,col),sy(m,m),wt(m,m)
 
c     ************
c
c     Subroutine subcg
c
c     This subroutine performs subspace minimization by means
c       of the conjugate gradient method. 
c
c     Given xcp, l, u, r, an index set ind that specifies the free
c       variables at xcp, and an l-BFGS matrix B (in terms of matrices
c       WY, WS, SY, WT, and scalars head, col, and theta) this subroutine
c       computes an approximate solution of the subspace problem
c
c     	min r'(Z'x-Z'xcp) + 1/2 (Z'x-Z'xcp)' (Z'B Z) (Z'x-Z'xcp)
c
c                     subject to l<=x<=u
c
c       by the conjugate gradient method. Only the free variables
c       specified by ind are changed in this subroutine.
c       
c     The subroutine statement is:
c
c       subroutine subcg(n,nsub,ind,l,u,nbd,x,r,p,q,m,ws,wy,sy,wt,
c                      theta,col,head,itercg,wv,iword,iprint)
c
c     where
c
c     n is an integer variable.
c       On entry n is total number of variables.
c       On exit n is unchanged.
c
c     nsub is an integer variable.
c       On entry nsub is the number of free variables.
c       On exit nsub is unchanged. 
c
c     ind is an integer array of dimension nsub.
c       On entry ind specifies the indices of the free variables.
c       On exit ind is unchanged. 
c
c     l is a double precision array of dimension n.
c       On entry l is the lower bound of x.
c       On exit l is unchanged.
c
c     u is a double precision array of dimension n.
c       On entry u is the upper bound of x.
c       On exit u is unchanged.
c
c     nbd is a integer array of dimension n.
c       On entry nbd represents the type of bounds imposed on the
c         variables:
c         nbd(i)=0 if x(i) is unbounded,
c                1 if x(i) has only a lower bound,
c                2 if x(i) has both lower and upper bounds, 
c                3 if x(i) has only upper bound.
c       On exit nbd is unchanged.
c
c     x is a double precision array of dimension n.
c       On entry x specifies the Generalized Cauchy point xcp.
c       On exit x is the subspace minimizer; i.e. x(i) = xcp(i)  
c          for i not in ind, and x(i) is determined by the
c          subspace minimization for i in ind. 
c
c     r is a double precision working array of dimension n.
c       On entry r is minus the reduced gradient of the model
c       at xcp.
c
c     p is a double precision working array of dimension nsub.
c
c     q is a double precision working array of dimension nsub.
c
c     m is an integer variable.
c       On entry m is the maximum number of variable metric 
c          corrections allowed in the limited memory matrix.
c       On exit m is unchanged.
c
c     ws, wy, sy, and wt are double precision arrays.
c       On entry they store the following information defining
c          the limited memory BFGS matrix:
c          ws, of dimension n x m, stores S, the matrix of s-vectors;
c          wy, of dimension n x m, stores Y, the matrix of y-vectors;
c          sy, of dimension m x m, stores S'Y;
c          wt, of dimension m x m, stores the Cholesky factorization
c                                  of (theta*S'S+LD^(-1)L'); see eq.
c                                  (2.26) in [1].
c       On exit ws, wy, sy, and wt are unchanged.
c
c     theta is a double precision variable;
c       On entry theta is the scaling factor in B_0 = theta I;
c       On exit theta is unchanged.
c
c     col is an integer variable;
c       On entry col is the number of variable metric corrections 
c          stored so far.
c       On exit col is unchanged.
c
c     head is an integer variable.
c       On entry head is the location of the 1st s-vector (or y-vector)
c          in S (or Y).
c       On exit head is unchanged.
c
c     itercg is an integer variable.
c       On entry itercg is unspecified.
c       On exit itercg specifies the number of cg iterations performed
c          during the subspace minimization.
c
c     wv is a double precision working array of dimension 2m.
c
c     iword is an integer variable.
c       On entry iword is unspecified.
c       On exit iword specifies the status of the subspace solution.
c  	  iword = 0 if the CG stopping rule is satisfied,
c                 1 if a bound is encountered.
c
c     iprint is an INTEGER variable that must be set by the user.
c       It controls the frequency and type of output generated:
c        iprint<0    no output is generated;
c        iprint=0    print only one line at the last iteration;
c        0<iprint<99 print also f and |proj g| every iprint iterations;
c        iprint=99   print details of every iteration except n-vectors;
c        iprint=100  print also the changes of active set and final x;
c        iprint>100  print details of every iteration including x and g;
c       When iprint > 0, the file iterate.dat will be created to
c                        summarize the iteration.
c
c     info is an integer variable.
c       On entry info is 0.
c       On exit info = 0       for normal return,
c                    = nonzero for abnormal return; this occurs when
c                              the subroutine bmv, which is called by
c                              subroutine subbv, detects a singular
c                              matrix.
c
c     References:
c     [1] R. Byrd, J. Nocedal and R. Schnabel "Representations of
c     Quasi-Newton Matrices and their use in Limited Memory Methods,
c     (1992), Northwestern University,  EECS Dept., Rep. NAM 06, to
c     appear in Mathematical Programming.
c
c     Subprograms called:
c
c       L-BFGS-B Library ... subbv.
c
c       Linpack ... dcopy, ddot.
c
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************

      double precision one,zero
      parameter(one=1.0d0,zero=0.0d0)

      integer i,k,inkeep,m2
      double precision ddot,rho1,rho2,beta,alpha1,alpha2,epsil,
     +                 pq,temp1,temp2,pi,pk
 
      if (nsub .le. 0) return
      itercg = 0
      if (iprint .ge. 99) write (6,1000)
 
c     Initialize CG direction p
 
      call dcopy(nsub,r,1,p,1)
 
c     Initialize rho2=r'r.
 
      rho2 = ddot(nsub,r,1,r,1)
 
c     Set tolerance.
 
      epsil = (min(1.0d-2,sqrt(rho2)))*rho2
      m2 = 2*m
 
c     CG iteration.
 
 111  if (rho2 .le. epsil .or. itercg .gt. m2) then
	 iword = 0
	 goto 222
      endif	
      itercg = itercg + 1
 
c     Compute q = (Z'B Z) p.
 
      if (col .gt. 0) then
         call subbv(n,m,nsub,ind,ws,wy,sy,wt,theta,col,head,p,q,wv,info)
	 if (info .ne. 0) return
      else
 	 do 10 i = 1, nsub
	    q(i) = theta*p(i)
  10     continue
      endif
      pq = ddot(nsub,p,1,q,1)
 
c     alpha2 = stepsize of exact minimizer along p.
 
      alpha2 = rho2/pq
 
c     alpha1 = largest feasible stepsize if smaller than alpha2.
 
      alpha1 = alpha2
      temp1 = alpha1		
      do 20 i = 1, nsub
	 k = ind(i)
         pi = p(i)
	 if (nbd(k) .ne. 0) then
   	    if (pi .lt. zero .and. nbd(k) .le. 2) then
	       temp2 = l(k) - x(k)
	       if (temp2 .ge. zero) then
		  temp1 = zero
	       else if (pi*alpha2 .lt. temp2) then
		  temp1 = temp2/pi
 	       endif
   	    else if (pi .gt. zero .and. nbd(k) .ge. 2) then
	       temp2 = u(k) - x(k)
	       if (temp2 .le. zero) then
		  temp1 = zero
	       else if (pi*alpha2 .gt. temp2) then
		  temp1 = temp2/pi
 	       endif
            endif
            if (temp1 .lt. alpha1) then
	       alpha1 = temp1
	       inkeep = i
            endif
         endif
  20  continue
 
c     The stepsize is the minimum of alpha1 and alpha2. 
 
      if (alpha2 .gt. alpha1) then
      	 pk = p(inkeep)
      	 k = ind(inkeep)
      	 if (pk .gt. zero) then
            x(k) = u(k)
            p(inkeep) = zero
     	 else if (pk .lt. zero) then
            x(k) = l(k)
	    p(inkeep) = zero
     	 endif
	 do 30 i = 1, nsub
	    k = ind(i)
	    x(k) = x(k) + alpha1*p(i)
  30  	 continue
	 iword=1
	 go to 222
      else
	 do 40 i = 1, nsub
            k = ind(i)
	    x(k) = x(k) + alpha2*p(i)
	    r(i) = r(i) - alpha2*q(i)
  40  	 continue
	 rho1 = rho2
	 rho2 = ddot(nsub,r,1,r,1)
	 beta = rho2/rho1
         do 50 i = 1, nsub
	    p(i) = r(i) + beta*p(i)
  50  	 continue
	 goto 111
      endif
 222  continue
      if (iprint .ge. 99) then
         write (6,1001) itercg
         if (iword .eq. 1) then
	    write (6,*) 'hit the boundary'
	 else
	    write (6,*) 'CG solution inside the box'
	 endif
      endif
      if(iprint .gt. 100) write (6,1002) (x(i),i=1,n)
      if (iprint .ge. 99) write (6,1003)
 1000 format (/,'----------------CG entered-----------------',/)
 1001 format ('CG iteration= ',i3)
 1002 format ('Subspace solution X =  ',/,(4x,1p,6(1x,d11.4)))
 1003 format (/,'----------------exit CG---------------------',/)

      return

      end

c====================== The end of subcg ===============================
 
      subroutine subdl(n,m,nact,ind,l,u,nbd,x,xcp,g,d,tlmda,
     +                   ws,wy,sy,yy,theta,col,head,sg,yg,
     +	                 iword,wv,wv2,wn,iprint,info)
 
      integer n,m,nact,col,head,iword,iprint,info,
     +        ind(n),nbd(n)
      double precision theta,
     +                 l(n),u(n),x(n),xcp(n),g(n),d(n),tlmda(n),
     +                 ws(n,m),wy(n,m),sy(m,m),yy(m,m),sg(m),yg(m),
     +                 wv(2*m),wv2(2*m),wn(2*m,2*m)

c     ************
c
c     Subroutine subdl
c
c     Given x, xcp, l, u, g, an index set that specifies
c	the active set at xcp, and an L-BFGS matrix H 
c	(in terms of WY, WS, SY, YY, head, col, and theta), 
c	this subroutine computes an approximate solution 
c	of the subspace problem
c
c     	(P)   min Q(z) = g'(z-x) + 1/2 (z-x)' B (z-x)  
c          
c             subject to   A'(z-x)=A'(xcp-x)          
c			   l<=z<=u
c
c       by a dual approach.  The solution is returned in xcp.
c       
c     The subroutine statement is:
c
c       subroutine subdl(n,m,nact,ind,l,u,nbd,x,xcp,g,d,tlmda,
c                          ws,wy,sy,yy,theta,col,head,sg,yg,
c     	                   iword,wv,wv2,wn,iprint,info)
c
c       where
c
c     n is an integer variable.
c       On entry n is the dimension of the problem.
c       On exit n is unchanged.
c
c     m is an integer variable.
c       On entry m is the maximum number of variable metric corrections
c         used to define the limited memory matrix.
c       On exit m is unchanged.
c
c     nact is an integer variable.
c       On entry nact is the number of active variables.
c       On exit nact is unchanged.
c
c     ind is an integer array of diminsion nact.
c       On entry ind specifies the indices of active variables.
c       On exit ind is unchanged.
c
c     l is a double precision array of dimension n.
c       On entry l is the lower bound of x.
c       On exit l is unchanged.
c
c     u is a double precision array of dimension n.
c       On entry u is the upper bound of x.
c       On exit u is unchanged.
c
c     nbd is a integer array of dimension n.
c       On entry nbd represents the type of bounds imposed on the
c         variables, and must be specified as follows:
c         nbd(i)=0 if x(i) is unbounded,
c                1 if x(i) has only lower bound,
c                2 if x(i) has both lower and upper bounds, and
c                3 if x(i) has only upper bound.
c       On exit nbd is unchanged.
c
c     x is a double precision array of dimension n.
c       On entry x specifies the current iterate x.  
c       On exit x is unchanged.
c
c     xcp is a double precision array of dimension n.
c       On entry xcp specifies the Cauchy point.  On exit xcp is an
c       approximate solution of (P).
c
c     g is a double precision array of dimension n.
c       On entry g stores the gradient of f at x.
c       On exit g is unchanged.
c
c     d is a double precision working array of dimension n.
c
c     tlmda is a double precision working array of dimension n. 
c
c     ws, wy, sy, and yy are double precision arrays;
c     theta is a double precision variable;
c     col is an integer variable;
c     head is an integer variable.
c       On entry they store the information defining the
c                                          limited memory BFGS matrix:
c         ws(n,m) stores S, a set of s-vectors;
c         wy(n,m) stores Y, a set of y-vectors;
c         sy(m,m) stores S'Y;
c         yy(m,m) stores Y'Y;
c         theta is the scaling factor specifying H_0 = I/theta;
c         col is the number of variable metric corrections stored;
c         head is the location of the 1st s- (or y-) vector in S (or Y).
c       On exit they are unchanged.
c
c     sg is a double precision array of dimension m.
c       On entry sg specifies S'g.
c       On exit sg is unchanged.
c
c     yg is a double precision array of dimension m.
c       On entry yg specifies Y'g.
c       On exit yg is unchanged.
c
c     iword is an integer variable.
c       On entry iword is unspecified.
c       On exit iword specifies the status of the subspace solution.
c	  iword = 0 if the solution is in the box,
c                 1 if some bound is encountered.
c
c     wv and  wv2 are double precision working arrays of dimension 2m.
c
c     wn is a double precision array of dimension 2m x 2m.
c       On entry the upper trangle of wn stores the LEL^T factorization
c         of the indefinite matrix
c              K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
c                  [L_a -R_z           theta*S'AA'S ]
c                                                    where E = [-I  0]
c                                                              [ 0  I]
c     Note that this procedure for computing d differs 
c     from the dual method described in the paper "A limited memory 
c     method for bound constrained optimization". 
c     The matrix K is equal to the matrix Mbar^[-1] Nbar described in 
c     secttion 5.3 that paper.
c       On exit wn is unchanged.
c
c     iprint is an INTEGER variable that must be set by the user.
c       It controls the frequency and type of output generated:
c        iprint<0    no output is generated;
c        iprint=0    print only one line at the last iteration;
c        0<iprint<99 print also f and |proj g| every iprint iterations;
c        iprint=99   print details of every iteration except n-vectors;
c        iprint=100  print also the changes of active set and final x;
c        iprint>100  print details of every iteration including x and g;
c       When iprint > 0, the file iterate.dat will be created to
c                        summarize the iteration.
c
c     info is an integer variable.
c       On entry info is unspecified.
c       On exit info = 0       for normal return,
c                    = nonzero for abnormal return
c                                  when the matrix K is ill-conditioned.
c
c     Subprograms called:
c
c
c       L-BFGS-B Library ... hmv.
c
c       Linpack ... dcopy, dtrsl, daxpy.
c
c      (In this version, the dual method is implemented by using g_k
c      instead of r^c, and W'g is formed by using S'g, Y'g in 'sblbfgs'.
c
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************

      double precision one,zero
      parameter(one=1.0d0,zero=0.0d0)

      integer i,j,k,ldt,col2,ibd,pointr,ihd
      double precision t1,t2,sum1,sum2,alpha,theta2,a1,a2,dk
 
      if (nact .ge. n) return
      if (iprint .ge. 99) write (6,1001)
 
c     wv=W'g=[Y'g/theta, S'g].
 
      do 10 i = 1, col
         wv(i) = yg(i)/theta
  10  continue
      call dcopy(col,sg,1,wv(col+1),1)
 
c     wv2=M(wv)=MW'g.
 
      call hmv(m,sy,yy,theta,col,wv,wv2,info)
      if (info .ne. 0) return 
 
      ihd = n - nact
      if(nact .gt. 0) then
c                          some variables are active.

c        Compute the rhs of the linear system for the  
c                                           Lagrange multiplier tlmda.
c          d=-(A'Hg+A'(x^c-x))=-A'g/theta-A'(x^c-x)-A'W(wv2), 
c            where A' is the "absolute" value of the gradient of
c                                               the active constraints.

         do 15 i = 1, nact
            k = ind(i + ihd)
            d(i) = -g(k)/theta - xcp(k) + x(k)
  15     continue

c          compute d:=d-A'[Y/theta S]*wv2 by a column-oriented approach.
         pointr=head 
         do 25 j=1, col
            a1 = wv2(j)/theta
            a2 = wv2(col + j)
            do 20 i = 1, nact
	       k = ind(i + ihd)
	       d(i) = d(i) - wy(k,pointr)*a1 - ws(k,pointr)*a2
  20   	    continue
            pointr = mod(pointr,m) + 1
  25     continue

c        Solve the linear system for tlmda using the
c                                       Sherman Morrison formula.
c          tlmda=(A'HA)^(-1)*d
c	        =theta*d-theta**2*A'WK^(-1)W'Ad.

c          first compute wv=W'Ad=[ Y'A*d/theta ]
c                                [ S'A*d       ].
         pointr = head 
         do 35 i = 1, col
            sum1 = zero
	    sum2 = zero
	    do 30 j = 1, nact
	       k = ind(j + ihd)
	       sum1 = sum1 + wy(k,pointr)*d(j)
	       sum2 = sum2 + ws(k,pointr)*d(j)
  30        continue
	    wv(i) = sum1/theta
	    wv(col + i) = sum2
            pointr = mod(pointr,m) + 1
  35     continue
 
c          compute wv:=K^(-1)wv.
 
         ldt = 2*m
         col2 = 2*col
         call dtrsl(wn,ldt,col2,wv,11,info)
         if (info .ne. 0) return 
	 do 37 i = 1, col
	     wv(i) = -wv(i)
  37     continue
         call dtrsl(wn,ldt,col2,wv,01,info)
         if (info .ne. 0) return 
 
c          tlmda = theta d - theta^2 A'W wv
         do 40 i = 1, nact
	    tlmda(i) = theta*d(i)
  40     continue
         theta2 = theta*theta
         pointr = head
         do 50 j = 1, col
	    a1 = wv(j)*theta
	    a2 = wv(j + col)*theta2
	    do 45 i = 1, nact
	       k = ind(i + ihd)
	       tlmda(i) =
     +                 tlmda(i) - wy(k,pointr)*a1 - ws(k,pointr)*a2
  45        continue
            pointr = mod(pointr,m) + 1
  50     continue
 
c        Solve the linear system of d.
c
c          d = -H(A tlmda+g)=-(A tlmda+g)/theta-W(MW'(A tlmda)+wv2) 
c                                                     [note wv2=MW'g].
 
c          wv=W'A tlmda.    
         pointr=head 
         do 60 j = 1, col
            wv(j) = zero
            wv(col + j) = zero
            do 55 i = 1, nact
	       k = ind(i + ihd)
	       wv(j) = wv(j) + wy(k,pointr)*tlmda(i)
	       wv(col + j) = wv(col+j) + ws(k,pointr)*tlmda(i)
  55   	    continue
            wv(j) = wv(j)/theta
            pointr = mod(pointr,m) + 1
  60     continue
 
c          wv=Mwv.
         call hmv(m,sy,yy,theta,col,wv,wv,info)
         if (info .ne. 0) return
 
c          wv2=wv+wv2.
         do 65 i = 1, 2*col
            wv2(i) = wv(i) + wv2(i)
  65     continue
      endif 
 
c       d:=-g/theta.
      do 70 i = 1, n
         d(i) = -g(i)/theta
  70  continue
      if(nact .gt. 0) then
c                          d:=d-(A tlmda)/theta.
         do 75 i = 1, nact
            k = ind(i + ihd)
            d(k) = d(k) - tlmda(i)/theta
  75     continue
      endif
    
c       d:=d-W(wv2).
      pointr=head 
      do 80 j = 1, col
         a1 = wv2(j)/theta
         a2 = wv2(col+j)
	 call daxpy(n,-a1,wy(1,pointr),1,d,1)
	 call daxpy(n,-a2,ws(1,pointr),1,d,1)
         pointr = mod(pointr,m) + 1
  80  continue
 
c     Backtrack to the feasible region.
 
c       new point tlmda=x+d.
      do 85 i = 1, n 
         tlmda(i) = x(i) + d(i)
         d(i) = tlmda(i) - xcp(i)
  85  continue
 
c       in case the floating point number d(i) .ne. zero
c                                                      for i \in A(x^c).
      do 90 i = 1, nact
	 k = ind(i + ihd)
         d(k) = zero
	 tlmda(k) = xcp(k)
  90  continue 
 
      alpha=one
      t1 = alpha   	
      do 95 k = 1, n
         dk = d(k) 
         if (nbd(k) .ne. 0) then
   	    if (dk .lt. zero .and. nbd(k) .le. 2) then
	       t2 = l(k) - xcp(k)
	       if (t2 .ge. zero) then
		  t1 = zero
	       else if (dk*alpha .lt. t2) then
		  t1 = t2/dk
 	       endif
   	    else if (dk .gt. zero .and. nbd(k) .ge. 2) then
	       t2 = u(k) - xcp(k)
	       if (t2 .le. zero) then
		  t1 = zero
	       else if (dk*alpha .gt. t2) then
		  t1 = t2/dk
 	       endif
            endif
            if (t1 .lt. alpha) then
	       alpha = t1
	       ibd = k
            endif
	 endif
  95  continue
 
c     The new point is a linear combination of tlmda and x^c
c       if tlmda is outside of the box; otherwise it is tlmda.
 
      if (alpha .lt. one) then
         call daxpy(n,alpha,d,1,xcp,1)
         if (d(ibd) .gt. zero) then
	    xcp(ibd) = u(ibd)
         else
	    xcp(ibd) = l(ibd)
         endif
         iword = 1
      else
         call dcopy(n,tlmda,1,xcp,1)
         iword = 0
      endif
 
      if (iprint .ge. 99) then
	 if (alpha .lt. one) then
            write (6,1101) alpha
         else
            write (6,*) 'solution inside the box'
	 endif	
	 if (iprint .gt.100) write (6,1201) (xcp(i),i=1,n)
      endif
      if (iprint .ge. 99) write (6,1301)
 1001 format (/,'----------------SUBDUAL entered-----------------',/)
 1101 format ( 'ALPHA = ',f7.5,' backtrack to the BOX')	
 1201 format ('Subspace solution X =  ',/,(4x,1p,6(1x,d11.4)))
 1301 format (/,'----------------exit SUBDUAL --------------------',/)
 
      return
 
      end

c====================== The end of subdl ===============================

      subroutine subsm(n,m,nsub,ind,l,u,nbd,x,d,ws,wy,theta,
     +                 col,head,iword,wv,wn,iprint,info)
 
      integer n,m,nsub,col,head,iword,iprint,info,
     +        ind(nsub),nbd(n)
      double precision theta,
     +                 l(n),u(n),x(n),d(n),
     +                 ws(n,m),wy(n,m),
     +                 wv(2*m),wn(2*m,2*m)

c     ************
c
c     Subroutine subsm
c
c     Given xcp, l, u, r, an index set that specifies
c	the active set at xcp, and an l-BFGS matrix B 
c	(in terms of WY, WS, SY, WT, head, col, and theta), 
c	this subroutine computes an approximate solution
c	of the subspace problem
c
c     	(P)   min Q(x) = r'(x-xcp) + 1/2 (x-xcp)' B (x-xcp)
c
c             subject to l<=x<=u
c	  	        x_i=xcp_i for all i in A(xcp)
c	              
c	along the subspace unconstrained Newton direction 
c	
c	   d = -(Z'BZ)^(-1) r.
c
c	The formula for the Newton direction, by the L-BFGS matrix and
c	the Sherman-Morrison formula, is
c
c	   d = (1/theta)r + (1/theta*2) Z'WK^(-1)W'Z r.
c 
c       where
c                 K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
c                     [L_a -R_z           theta*S'AA'S ]
c     Note that this procedure for computing d differs 
c     from that described in the paper "A limited memory method
c     for bound constrained optimization". The matrix K is equal to the 
c     matrix M^[-1]N in that paper.
c
c     The subroutine statement is:
c
c       subroutine subsm(n,m,nsub,ind,l,u,nbd,x,d,ws,wy,theta,
c                        col,head,iword,wv,wn,iprint,info)
c
c       where
c
c     n is an integer variable.
c       On entry n is the dimension of the problem.
c       On exit n is unchanged.
c
c     m is an integer variable.
c       On entry m is the maximum number of variable metric corrections
c         used to define the limited memory matrix.
c       On exit m is unchanged.
c
c     nsub is an integer variable.
c       On entry nsub is the number of free variables.
c       On exit nsub is unchanged.
c
c     ind is an integer array of dimension nsub.
c       On entry ind specifies the coordinate indices of free variables.
c       On exit ind is unchanged.
c
c     l is a double precision array of dimension n.
c       On entry l is the lower bound of x.
c       On exit l is unchanged.
c
c     u is a double precision array of dimension n.
c       On entry u is the upper bound of x.
c       On exit u is unchanged.
c
c     nbd is a integer array of dimension n.
c       On entry nbd represents the type of bounds imposed on the
c         variables, and must be specified as follows:
c         nbd(i)=0 if x(i) is unbounded,
c                1 if x(i) has only lower bound,
c                2 if x(i) has both lower and upper bounds, and
c                3 if x(i) has only upper bound.
c       On exit nbd is unchanged.
c
c     x is a double precision array of dimension n.
c       On entry x specifies the Cauchy point xcp. 
c       On exit x(i) is the minimizer of Q over the subspace of
c                                                        free variables. 
c
c     d is a double precision array of dimension n.
c       On entry d is the reduced gradient of Q at xcp.
c       On exit d is the Newton direction of Q. 
c
c     ws and wy are double precision arrays;
c     theta is a double precision variable;
c     col is an integer variable;
c     head is an integer variable.
c       On entry they store the information defining the
c                                          limited memory BFGS matrix:
c         ws(n,m) stores S, a set of s-vectors;
c         wy(n,m) stores Y, a set of y-vectors;
c         theta is the scaling factor specifying B_0 = theta I;
c         col is the number of variable metric corrections stored;
c         head is the location of the 1st s- (or y-) vector in S (or Y).
c       On exit they are unchanged.
c
c     iword is an integer variable.
c       On entry iword is unspecified.
c       On exit iword specifies the status of the subspace solution.
c         iword = 0 if the solution is in the box,
c                 1 if some bound is encountered.
c
c     wv are double precision working arrays of dimension 2m.
c
c     wn is a double precision array of dimension 2m x 2m.
c       On entry the upper triangle of wn stores the LEL^T factorization
c         of the indefinite matrix
c              K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
c                  [L_a -R_z           theta*S'AA'S ]
c                                                    where E = [-I  0]
c                                                              [ 0  I]
c       On exit wn is unchanged.
c
c     iprint is an INTEGER variable that must be set by the user.
c       It controls the frequency and type of output generated:
c        iprint<0    no output is generated;
c        iprint=0    print only one line at the last iteration;
c        0<iprint<99 print also f and |proj g| every iprint iterations;
c        iprint=99   print details of every iteration except n-vectors;
c        iprint=100  print also the changes of active set and final x;
c        iprint>100  print details of every iteration including x and g;
c       When iprint > 0, the file iterate.dat will be created to
c                        summarize the iteration.
c
c     info is an integer variable.
c       On entry info is unspecified.
c       On exit info = 0       for normal return,
c                    = nonzero for abnormal return 
c                                  when the matrix K is ill-conditioned.
c
c     Subprograms called:
c
c       Linpack dtrsl.
c
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************

      double precision one,zero
      parameter(one=1.0d0,zero=0.0d0)

      integer pointr,m2,col2,ibd,jy,js,i,j,k
      double precision alpha,dk,temp1,temp2
 
      if (nsub .le. 0) return
      if (iprint .ge. 99) write (6,1001)

c     Compute wv = W'Zd.

      pointr = head 
      do 20 i = 1, col
     	 temp1 = zero
	 temp2 = zero
	 do 10 j = 1, nsub
	    k = ind(j)
	    temp1 = temp1 + wy(k,pointr)*d(j)
	    temp2 = temp2 + ws(k,pointr)*d(j)
  10     continue
	 wv(i) = temp1
	 wv(col + i) = theta*temp2
	 pointr = mod(pointr,m) + 1
  20  continue
 
c     Compute wv:=K^(-1)wv.

      m2 = 2*m
      col2 = 2*col
      call dtrsl(wn,m2,col2,wv,11,info)
      if (info .ne. 0) return
      do 25 i = 1, col
	 wv(i) = -wv(i)
  25     continue
      call dtrsl(wn,m2,col2,wv,01,info)
      if (info .ne. 0) return
 
c     Compute d = (1/theta)d + (1/theta**2)Z'W wv.
 
      pointr = head
      do 40 jy = 1, col
         js = col + jy
	 do 30 i = 1, nsub
	    k = ind(i)
	    d(i) = d(i) + wy(k,pointr)*wv(jy)/theta     
     +                  + ws(k,pointr)*wv(js)
  30     continue
	 pointr = mod(pointr,m) + 1
  40  continue
      do 50 i = 1, nsub
	 d(i) = d(i)/theta
  50  continue
 
c     Backtrack to the feasible region.
 
      alpha = one
      temp1 = alpha	
      do 60 i = 1, nsub
	 k = ind(i)
         dk = d(i)
	 if (nbd(k) .ne. 0) then
   	    if (dk .lt. zero .and. nbd(k) .le. 2) then
	       temp2 = l(k) - x(k)
	       if (temp2 .ge. zero) then
		  temp1 = zero
	       else if (dk*alpha .lt. temp2) then
		  temp1 = temp2/dk
 	       endif
   	    else if (dk .gt. zero .and. nbd(k) .ge. 2) then
	       temp2 = u(k) - x(k)
	       if (temp2 .le. zero) then
		  temp1 = zero
	       else if (dk*alpha .gt. temp2) then
		  temp1 = temp2/dk
 	       endif
            endif
            if (temp1 .lt. alpha) then
	       alpha = temp1
	       ibd = i
            endif
         endif
  60  continue
 
      if (alpha .lt. one) then
      	 dk = d(ibd)
      	 k = ind(ibd)
      	 if (dk .gt. zero) then
            x(k) = u(k)
            d(ibd) = zero
     	 else if (dk .lt. zero) then
            x(k) = l(k)
	    d(ibd) = zero
     	 endif
      endif
      do 70 i = 1, nsub
	 k = ind(i)
	 x(k) = x(k) + alpha*d(i)
  70  continue
 
      if (iprint .ge. 99) then
	 if (alpha .lt. one) then
            write (6,1002) alpha
         else
            write (6,*) 'SM solution inside the box'
	 end if	
	 if (iprint .gt.100) write (6,1003) (x(i),i=1,n)
      endif
 
      if (alpha .lt. one) then
         iword = 1
      else
         iword = 0
      endif 
      if (iprint .ge. 99) write (6,1004)

 1001 format (/,'----------------SUBSM entered-----------------',/)
 1002 format ( 'ALPHA = ',f7.5,' backtrack to the BOX')	
 1003 format ('Subspace solution X =  ',/,(4x,1p,6(1x,d11.4)))
 1004 format (/,'----------------exit SUBSM --------------------',/)

      return

      end
      
c====================== The end of subsm ===============================

      subroutine sytimg(n,m,g,ws,wy,sg,sgo,yg,ygo,col,head,iupdat)
 
      integer n,m,col,head,iupdat
      double precision g(n),ws(n,m),wy(n,m),sg(m),yg(m),sgo(m),ygo(m)

c     ************
c
c     Subroutine sytimg
c
c       This subroutine computes S'g and Y'g.
c
c     Subprograms called:
c
c       Linpack ... dcopy, ddot.
c
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c
c     ************

      integer i,pointr
      double precision ddot
     
      i = 1
      if (iupdat .gt. m) i = 2
      call dcopy(col-1, sg(i),1,sgo,1)
      call dcopy(col-1, yg(i),1,ygo,1)
      pointr = head
      do 56 i = 1, col
         sg(i) = ddot(n,ws(1,pointr),1,g,1)
         yg(i) = ddot(n,wy(1,pointr),1,g,1)
         pointr = mod(pointr,m) + 1
  56  continue
 
      return

      end

c======================= The end of sytimg =============================

      subroutine heapst(n,t,iorder,iheap)
      integer iheap,n,iorder(n)
      double precision t(n)
  
c     ************
c
c     Subroutine heapst 
c
c     This subroutine sorts out the least element of t, and puts the
c       remaining elements of t in a heap.
c 
c     The subroutine statement is:
c
c       subroutine heapst(n,t,iorder,iheap)
c
c     where
c
c     n is an integer variable.
c       On entry n is the dimension of the arrays t and iorder.
c       On exit n is unchanged.
c
c     t is an double precision array of dimension n.
c       On entry t stores the elements to be sorted,
c       On exit t(n) stores the least elements of t, and t(1) to t(n-1)
c         stores the remaining elements in the form of a heap.
c
c     iorder is an integer array of dimension n.
c       On entry iorder(i) is the index of t(i).
c       On exit iorder(i) still is the index of t(i), but iorder may be
c         permuted in accordance with t.
c
c     iheap is an integer variable specifing the task.
c       On entry iheap should be set as follows:
c         iheap .eq. 0 if t(1) to t(n) is not in the form of a heap,
c         iheap .ne. 0 if otherwise.
c       On exit iheap is unchanged.
c
c
c     References:
c       Algorithm 232 of CACM (J. W. J. Williams): HEAPSORT.
c
c                           *  *  *
c
c     NEOS, November 1994.
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c     ************
  
      integer i,j,k,indxin,indxou
      double precision ddum,out

      if (iheap .eq. 0) then

c        Rearrange the elements t(1) to t(n) to form a heap.

         do 20 k = 2, n
            ddum  = t(k)
            indxin = iorder(k)

c           Add ddum to the heap.
            i = k
   10       continue
            if (i.gt.1) then
               j = i/2
               if (ddum .lt. t(j)) then
                  t(i) = t(j)
                  iorder(i) = iorder(j)
                  i = j
                  goto 10 
               endif  
            endif  
            t(i) = ddum
            iorder(i) = indxin
   20    continue
      endif
 
c     Assign to 'out' the value of t(1), the least member of the heap,
c        and rearranges the remaining members to form a heap as
c        elements 1 to n-1 of t.
 
      if (n .gt. 1) then
         i = 1
         out = t(1)
         indxou = iorder(1)
         ddum  = t(n)
         indxin  = iorder(n)

c        Restore the heap 
   30    continue
         j = i+i
         if (j .le. n-1) then
            if (t(j+1) .lt. t(j)) j = j+1
            if (t(j) .lt. ddum ) then
               t(i) = t(j)
               iorder(i) = iorder(j)
               i = j
               goto 30
            endif 
         endif 
         t(i) = ddum
         iorder(i) = indxin
 
c     Put the least member in t(n). 

         t(n) = out
         iorder(n) = indxou
      endif 

      return

      end

c====================== The end of heapst ==============================

      subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax,
     +                  task,isave,dsave)
      character*(*) task
      integer isave(2)
      double precision f,g,stp,ftol,gtol,xtol,stpmin,stpmax
      double precision dsave(13)
c     **********
c
c     Subroutine dcsrch
c
c     This subroutine finds a step that satisfies a sufficient
c     decrease condition and a curvature condition.
c
c     Each call of the subroutine updates an interval with 
c     endpoints stx and sty. The interval is initially chosen 
c     so that it contains a minimizer of the modified function
c
c           psi(stp) = f(stp) - f(0) - ftol*stp*f'(0).
c
c     If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the
c     interval is chosen so that it contains a minimizer of f. 
c
c     The algorithm is designed to find a step that satisfies 
c     the sufficient decrease condition 
c
c           f(stp) <= f(0) + ftol*stp*f'(0),
c
c     and the curvature condition
c
c           abs(f'(stp)) <= gtol*abs(f'(0)).
c
c     If ftol is less than gtol and if, for example, the function
c     is bounded below, then there is always a step which satisfies
c     both conditions. 
c
c     If no step can be found that satisfies both conditions, then 
c     the algorithm stops with a warning. In this case stp only 
c     satisfies the sufficient decrease condition.
c
c     A typical invocation of dcsrch has the following outline:
c
c     task = 'START'
c  10 continue
c        call dcsrch( ... )
c        if (task .eq. 'FG') then
c           Evaluate the function and the gradient at stp 
c           goto 10
c           end if
c
c     NOTE: The user must no alter work arrays between calls.
c
c     The subroutine statement is
c
c        subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax,
c                          task,isave,dsave)
c     where
c
c       f is a double precision variable.
c         On initial entry f is the value of the function at 0.
c            On subsequent entries f is the value of the 
c            function at stp.
c         On exit f is the value of the function at stp.
c
c	g is a double precision variable.
c         On initial entry g is the derivative of the function at 0.
c            On subsequent entries g is the derivative of the 
c            function at stp.
c         On exit g is the derivative of the function at stp.
c
c	stp is a double precision variable. 
c         On entry stp is the current estimate of a satisfactory 
c            step. On initial entry, a positive initial estimate 
c            must be provided. 
c         On exit stp is the current estimate of a satisfactory step
c            if task = 'FG'. If task = 'CONV' then stp satisfies
c            the sufficient decrease and curvature condition.
c
c       ftol is a double precision variable.
c         On entry ftol specifies a nonnegative tolerance for the 
c            sufficient decrease condition.
c         On exit ftol is unchanged.
c
c       gtol is a double precision variable.
c         On entry gtol specifies a nonnegative tolerance for the 
c            curvature condition. 
c         On exit gtol is unchanged.
c
c	xtol is a double precision variable.
c         On entry xtol specifies a nonnegative relative tolerance
c            for an acceptable step. The subroutine exits with a
c            warning if the relative difference between sty and stx
c            is less than xtol.
c         On exit xtol is unchanged.
c
c	stpmin is a double precision variable.
c         On entry stpmin is a nonnegative lower bound for the step.
c         On exit stpmin is unchanged.
c
c	stpmax is a double precision variable.
c         On entry stpmax is a nonnegative upper bound for the step.
c         On exit stpmax is unchanged.
c
c       task is a character variable of length at least 60.
c         On initial entry task must be set to 'START'.
c         On exit task indicates the required action:
c
c            If task(1:2) = 'FG' then evaluate the function and 
c            derivative at stp and call dcsrch again.
c
c            If task(1:4) = 'CONV' then the search is successful.
c
c            If task(1:4) = 'WARN' then the subroutine is not able
c            to satisfy the convergence conditions. The exit value of
c            stp contains the best point found during the search.
c
c            If task(1:5) = 'ERROR' then there is an error in the
c            input arguments.
c
c         On exit with convergence, a warning or an error, the
c            variable task contains additional information.
c
c       isave is an integer work array of dimension 2.
c         
c       dsave is a double precision work array of dimension 13.
c
c     Subprograms called
c
c	MINPACK-2 ... dcstep
c
c     MINPACK-1 Project. June 1983.
c     Argonne National Laboratory. 
c     Jorge J. More' and David J. Thuente.
c
c     MINPACK-2 Project. October 1993.
c     Argonne National Laboratory and University of Minnesota. 
c     Brett M. Averick, Richard G. Carter, and Jorge J. More'. 
c
c     **********
      double precision zero,p5,p66
      parameter(zero=0.0d0,p5=0.5d0,p66=0.66d0)
      double precision xtrapl,xtrapu
      parameter(xtrapl=1.1d0,xtrapu=4.0d0)

      logical brackt
      integer stage
      double precision finit,ftest,fm,fx,fxm,fy,fym,ginit,gtest,
     +       gm,gx,gxm,gy,gym,stx,sty,stmin,stmax,width,width1

c     Initialization block.

      if (task(1:5) .eq. 'START') then

c        Check the input arguments for errors.

         if (stp .lt. stpmin) task = 'ERROR: STP .LT. STPMIN'
         if (stp .gt. stpmax) task = 'ERROR: STP .GT. STPMAX'
         if (g .ge. zero) task = 'ERROR: INITIAL G .GE. ZERO'
         if (ftol .lt. zero) task = 'ERROR: FTOL .LT. ZERO'
         if (gtol .lt. zero) task = 'ERROR: GTOL .LT. ZERO'
         if (xtol .lt. zero) task = 'ERROR: XTOL .LT. ZERO'
         if (stpmin .lt. zero) task = 'ERROR: STPMIN .LT. ZERO'
         if (stpmax .lt. stpmin) task = 'ERROR: STPMAX .LT. STPMIN'

c        Exit if there are errors on input.

         if (task(1:5) .eq. 'ERROR') return

c        Initialize local variables.

         brackt = .false.
         stage = 1
         finit = f
         ginit = g
         gtest = ftol*ginit
         width = stpmax - stpmin
         width1 = width/p5

c        The variables stx, fx, gx contain the values of the step, 
c        function, and derivative at the best step. 
c        The variables sty, fy, gy contain the value of the step, 
c        function, and derivative at sty.
c        The variables stp, f, g contain the values of the step, 
c        function, and derivative at stp.

         stx = zero
         fx = finit
         gx = ginit
         sty = zero
         fy = finit
         gy = ginit
         stmin = zero
         stmax = stp + xtrapu*stp
         task = 'FG'

         goto 1000

      else

c        Restore local variables.

         if (isave(1) .eq. 1) then
            brackt = .true.
         else
            brackt = .false.
         endif
         stage = isave(2) 
         ginit = dsave(1) 
         gtest = dsave(2) 
         gx = dsave(3) 
         gy = dsave(4) 
         finit = dsave(5) 
         fx = dsave(6) 
         fy = dsave(7) 
         stx = dsave(8) 
         sty = dsave(9) 
         stmin = dsave(10) 
         stmax = dsave(11) 
         width = dsave(12) 
         width1 = dsave(13) 

      endif

c     If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the
c     algorithm enters the second stage.

      ftest = finit + stp*gtest
      if (stage .eq. 1 .and. f .le. ftest .and. g .ge. zero) 
     +   stage = 2

c     Test for warnings.

      if (brackt .and. (stp .le. stmin .or. stp .ge. stmax))
     +   task = 'WARNING: ROUNDING ERRORS PREVENT PROGRESS'
      if (brackt .and. stmax - stmin .le. xtol*stmax) 
     +   task = 'WARNING: XTOL TEST SATISFIED'
      if (stp .eq. stpmax .and. f .le. ftest .and. g .le. gtest) 
     +   task = 'WARNING: STP = STPMAX'
      if (stp .eq. stpmin .and. (f .gt. ftest .or. g .ge. gtest)) 
     +   task = 'WARNING: STP = STPMIN'

c     Test for convergence.

      if (f .le. ftest .and. abs(g) .le. gtol*(-ginit)) 
     +   task = 'CONVERGENCE'

c     Test for termination.

      if (task(1:4) .eq. 'WARN' .or. task(1:4) .eq. 'CONV') goto 1000

c     A modified function is used to predict the step during the
c     first stage if a lower function value has been obtained but 
c     the decrease is not sufficient.

      if (stage .eq. 1 .and. f .le. fx .and. f .gt. ftest) then

c        Define the modified function and derivative values.

         fm = f - stp*gtest
         fxm = fx - stx*gtest
         fym = fy - sty*gtest
         gm = g - gtest
         gxm = gx - gtest
         gym = gy - gtest

c        Call dcstep to update stx, sty, and to compute the new step.

         call dcstep(stx,fxm,gxm,sty,fym,gym,stp,fm,gm,
     +               brackt,stmin,stmax)

c        Reset the function and derivative values for f.

         fx = fxm + stx*gtest
         fy = fym + sty*gtest
         gx = gxm + gtest
         gy = gym + gtest

      else

c       Call dcstep to update stx, sty, and to compute the new step.

        call dcstep(stx,fx,gx,sty,fy,gy,stp,f,g,
     +              brackt,stmin,stmax)

      endif

c     Decide if a bisection step is needed.

      if (brackt) then
         if (abs(sty-stx) .ge. p66*width1) stp = stx + p5*(sty - stx)
         width1 = width
         width = abs(sty-stx)
      endif

c     Set the minimum and maximum steps allowed for stp.

      if (brackt) then
         stmin = min(stx,sty)
         stmax = max(stx,sty)
      else
         stmin = stp + xtrapl*(stp - stx)
         stmax = stp + xtrapu*(stp - stx)
      endif
 
c     Force the step to be within the bounds stpmax and stpmin.
 
      stp = max(stp,stpmin)
      stp = min(stp,stpmax)

c     If further progress is not possible, let stp be the best
c     point obtained during the search.

      if (brackt .and. (stp .le. stmin .or. stp .ge. stmax)
     +   .or. (brackt .and. stmax-stmin .le. xtol*stmax)) stp = stx

c     Obtain another function and derivative.

      task = 'FG'

 1000 continue

c     Save local variables.

      if (brackt) then
         isave(1) = 1
      else
         isave(1) = 0
      endif
      isave(2) = stage
      dsave(1) =  ginit
      dsave(2) =  gtest
      dsave(3) =  gx
      dsave(4) =  gy
      dsave(5) =  finit
      dsave(6) =  fx
      dsave(7) =  fy
      dsave(8) =  stx
      dsave(9) =  sty
      dsave(10) = stmin
      dsave(11) = stmax
      dsave(12) = width
      dsave(13) = width1

      end
      subroutine dcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt,
     +                  stpmin,stpmax)
      logical brackt
      double precision stx,fx,dx,sty,fy,dy,stp,fp,dp,stpmin,stpmax
c     **********
c
c     Subroutine dcstep
c
c     This subroutine computes a safeguarded step for a search
c     procedure and updates an interval that contains a step that
c     satisfies a sufficient decrease and a curvature condition.
c
c     The parameter stx contains the step with the least function
c     value. If brackt is set to .true. then a minimizer has
c     been bracketed in an interval with endpoints stx and sty.
c     The parameter stp contains the current step. 
c     The subroutine assumes that if brackt is set to .true. then
c
c           min(stx,sty) < stp < max(stx,sty),
c
c     and that the derivative at stx is negative in the direction 
c     of the step.
c
c     The subroutine statement is
c
c       subroutine dcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt,
c                         stpmin,stpmax)
c
c     where
c
c       stx is a double precision variable.
c         On entry stx is the best step obtained so far and is an
c            endpoint of the interval that contains the minimizer. 
c         On exit stx is the updated best step.
c
c       fx is a double precision variable.
c         On entry fx is the function at stx.
c         On exit fx is the function at stx.
c
c       dx is a double precision variable.
c         On entry dx is the derivative of the function at 
c            stx. The derivative must be negative in the direction of 
c            the step, that is, dx and stp - stx must have opposite 
c            signs.
c         On exit dx is the derivative of the function at stx.
c
c       sty is a double precision variable.
c         On entry sty is the second endpoint of the interval that 
c            contains the minimizer.
c         On exit sty is the updated endpoint of the interval that 
c            contains the minimizer.
c
c       fy is a double precision variable.
c         On entry fy is the function at sty.
c         On exit fy is the function at sty.
c
c       dy is a double precision variable.
c         On entry dy is the derivative of the function at sty.
c         On exit dy is the derivative of the function at the exit sty.
c
c       stp is a double precision variable.
c         On entry stp is the current step. If brackt is set to .true.
c            then on input stp must be between stx and sty. 
c         On exit stp is a new trial step.
c
c       fp is a double precision variable.
c         On entry fp is the function at stp
c         On exit fp is unchanged.
c
c       dp is a double precision variable.
c         On entry dp is the the derivative of the function at stp.
c         On exit dp is unchanged.
c
c       brackt is an logical variable.
c         On entry brackt specifies if a minimizer has been bracketed.
c            Initially brackt must be set to .false.
c         On exit brackt specifies if a minimizer has been bracketed.
c            When a minimizer is bracketed brackt is set to .true.
c
c       stpmin is a double precision variable.
c         On entry stpmin is a lower bound for the step.
c         On exit stpmin is unchanged.
c
c       stpmax is a double precision variable.
c         On entry stpmax is an upper bound for the step.
c         On exit stpmax is unchanged.
c
c     MINPACK-1 Project. June 1983
c     Argonne National Laboratory. 
c     Jorge J. More' and David J. Thuente.
c
c     MINPACK-2 Project. October 1993.
c     Argonne National Laboratory and University of Minnesota. 
c     Brett M. Averick and Jorge J. More'.
c
c     **********
      double precision zero,p66,two,three
      parameter(zero=0.0d0,p66=0.66d0,two=2.0d0,three=3.0d0)
   
      double precision gamma,p,q,r,s,sgnd,stpc,stpf,stpq,theta

      sgnd = dp*(dx/abs(dx))

c     First case: A higher function value. The minimum is bracketed. 
c     If the cubic step is closer to stx than the quadratic step, the 
c     cubic step is taken, otherwise the average of the cubic and 
c     quadratic steps is taken.

      if (fp .gt. fx) then
         theta = three*(fx - fp)/(stp - stx) + dx + dp
         s = max(abs(theta),abs(dx),abs(dp))
         gamma = s*sqrt((theta/s)**2 - (dx/s)*(dp/s))
         if (stp .lt. stx) gamma = -gamma
         p = (gamma - dx) + theta
         q = ((gamma - dx) + gamma) + dp
         r = p/q
         stpc = stx + r*(stp - stx)
         stpq = stx + ((dx/((fx - fp)/(stp - stx) + dx))/two)*
     +                                                       (stp - stx)
         if (abs(stpc-stx) .lt. abs(stpq-stx)) then
            stpf = stpc
         else
            stpf = stpc + (stpq - stpc)/two
         endif
         brackt = .true.

c     Second case: A lower function value and derivatives of opposite 
c     sign. The minimum is bracketed. If the cubic step is farther from 
c     stp than the secant step, the cubic step is taken, otherwise the 
c     secant step is taken.

      else if (sgnd .lt. zero) then
         theta = three*(fx - fp)/(stp - stx) + dx + dp
         s = max(abs(theta),abs(dx),abs(dp))
         gamma = s*sqrt((theta/s)**2 - (dx/s)*(dp/s))
         if (stp .gt. stx) gamma = -gamma
         p = (gamma - dp) + theta
         q = ((gamma - dp) + gamma) + dx
         r = p/q
         stpc = stp + r*(stx - stp)
         stpq = stp + (dp/(dp - dx))*(stx - stp)
         if (abs(stpc-stp) .gt. abs(stpq-stp)) then
            stpf = stpc
         else
            stpf = stpq
         endif
         brackt = .true.

c     Third case: A lower function value, derivatives of the same sign,
c     and the magnitude of the derivative decreases.

      else if (abs(dp) .lt. abs(dx)) then

c        The cubic step is computed only if the cubic tends to infinity 
c        in the direction of the step or if the minimum of the cubic
c        is beyond stp. Otherwise the cubic step is defined to be the 
c        secant step.

         theta = three*(fx - fp)/(stp - stx) + dx + dp
         s = max(abs(theta),abs(dx),abs(dp))

c        The case gamma = 0 only arises if the cubic does not tend
c        to infinity in the direction of the step.

         gamma = s*sqrt(max(zero,(theta/s)**2-(dx/s)*(dp/s)))
         if (stp .gt. stx) gamma = -gamma
         p = (gamma - dp) + theta
         q = (gamma + (dx - dp)) + gamma
         r = p/q
         if (r .lt. zero .and. gamma .ne. zero) then
            stpc = stp + r*(stx - stp)
         else if (stp .gt. stx) then
            stpc = stpmax
         else
            stpc = stpmin
         endif
         stpq = stp + (dp/(dp - dx))*(stx - stp)

         if (brackt) then

c           A minimizer has been bracketed. If the cubic step is 
c           closer to stp than the secant step, the cubic step is 
c           taken, otherwise the secant step is taken.

            if (abs(stpc-stp) .lt. abs(stpq-stp)) then
               stpf = stpc
            else
               stpf = stpq
            endif
            if (stp .gt. stx) then
               stpf = min(stp+p66*(sty-stp),stpf)
            else
               stpf = max(stp+p66*(sty-stp),stpf)
            endif
         else

c           A minimizer has not been bracketed. If the cubic step is 
c           farther from stp than the secant step, the cubic step is 
c           taken, otherwise the secant step is taken.

            if (abs(stpc-stp) .gt. abs(stpq-stp)) then
               stpf = stpc
            else
               stpf = stpq
            endif
            stpf = min(stpmax,stpf)
            stpf = max(stpmin,stpf)
         endif

c     Fourth case: A lower function value, derivatives of the same sign, 
c     and the magnitude of the derivative does not decrease. If the 
c     minimum is not bracketed, the step is either stpmin or stpmax, 
c     otherwise the cubic step is taken.

      else
         if (brackt) then
            theta = three*(fp - fy)/(sty - stp) + dy + dp
            s = max(abs(theta),abs(dy),abs(dp))
            gamma = s*sqrt((theta/s)**2 - (dy/s)*(dp/s))
            if (stp .gt. sty) gamma = -gamma
            p = (gamma - dp) + theta
            q = ((gamma - dp) + gamma) + dy
            r = p/q
            stpc = stp + r*(sty - stp)
            stpf = stpc
         else if (stp .gt. stx) then
            stpf = stpmax
         else
            stpf = stpmin
         endif
      endif

c     Update the interval which contains a minimizer.

      if (fp .gt. fx) then
         sty = stp
         fy = fp
         dy = dp
      else
         if (sgnd .lt. zero) then
            sty = stx
            fy = fx
            dy = dx
         endif
         stx = stp
         fx = fp
         dx = dp
      endif

c     Compute the new step.

      stp = stpf

      end
      subroutine timer(ttime)
      double precision ttime
c     *********
c
c     Subroutine timer
c
c     This subroutine is used to determine user time. In a typical 
c     application, the user time for a code segment requires calls 
c     to subroutine timer to determine the initial and final time.
c
c     The subroutine statement is
c
c       subroutine timer(ttime)
c
c     where
c
c       ttime is an output variable which specifies the user time.
c
c     Argonne National Laboratory and University of Minnesota.
c     MINPACK-2 Project.
c
c     Modified October 1990 by Brett M. Averick.
c
c     **********
      real temp
      real tarray(2)
      real etime

c     The first element of the array tarray specifies user time

      !temp = etime(tarray) 

      ttime = 0.0 !dble(tarray(1))
 
      return

      end
      double precision function dpmeps()
c     **********
c
c     Subroutine dpeps
c
c     This subroutine computes the machine precision parameter
c     dpmeps as the smallest floating point number such that
c     1 + dpmeps differs from 1.
c
c     This subroutine is based on the subroutine machar described in
c
c     W. J. Cody,
c     MACHAR: A subroutine to dynamically determine machine parameters,
c     ACM Transactions on Mathematical Software, 14, 1988, pages 303-311.
c
c     The subroutine statement is:
c
c       subroutine dpeps(dpmeps)
c
c     where
c
c       dpmeps is a double precision variable.
c         On entry dpmeps need not be specified.
c         On exit dpmeps is the machine precision.
c
c     MINPACK-2 Project. February 1991.
c     Argonne National Laboratory and University of Minnesota.
c     Brett M. Averick.
c
c     *******
      integer i,ibeta,irnd,it,itemp,negep
      double precision a,b,beta,betain,betah,temp,tempa,temp1,
     +       zero,one,two
      data zero,one,two /0.0d0,1.0d0,2.0d0/
 
c     determine ibeta, beta ala malcolm.

      a = one
      b = one
   10 continue
         a = a + a
         temp = a + one
         temp1 = temp - a
      if (temp1 - one .eq. zero) go to 10
   20 continue
         b = b + b
         temp = a + b
         itemp = int(temp - a)
      if (itemp .eq. 0) go to 20
      ibeta = itemp
      beta = dble(ibeta)

c     determine it, irnd.

      it = 0
      b = one
   30 continue
         it = it + 1
         b = b * beta
         temp = b + one
         temp1 = temp - b
      if (temp1 - one .eq. zero) go to 30
      irnd = 0
      betah = beta/two
      temp = a + betah
      if (temp - a .ne. zero) irnd = 1
      tempa = a + beta
      temp = tempa + betah
      if ((irnd .eq. 0) .and. (temp - tempa .ne. zero)) irnd = 2

c     determine dpmeps.

      negep = it + 3
      betain = one/beta
      a = one
      do 40 i = 1, negep
         a = a*betain
   40 continue
   50 continue
        temp = one + a
        if (temp - one .ne. zero) go to 60
        a = a*beta
        go to  50
   60 continue
      dpmeps = a
      if ((ibeta .eq. 2) .or. (irnd .eq. 0)) go to 70
      a = (a*(one + a))/two
      temp = one + a
      if (temp - one .ne. zero) dpmeps = a

   70 return

      end
      subroutine daxpy(n,da,dx,incx,dy,incy)
c
c     constant times a vector plus a vector.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(1),dy(1),da
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if (da .eq. 0.0d0) return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dy(iy) = dy(iy) + da*dx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,4)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dy(i) = dy(i) + da*dx(i)
   30 continue
      if( n .lt. 4 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,4
        dy(i) = dy(i) + da*dx(i)
        dy(i + 1) = dy(i + 1) + da*dx(i + 1)
        dy(i + 2) = dy(i + 2) + da*dx(i + 2)
        dy(i + 3) = dy(i + 3) + da*dx(i + 3)
   50 continue
      return
      end
      subroutine dcopy(n,dx,incx,dy,incy)
c
c     copies a vector, x, to a vector, y.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(1),dy(1)
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dy(iy) = dx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,7)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dy(i) = dx(i)
   30 continue
      if( n .lt. 7 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,7
        dy(i) = dx(i)
        dy(i + 1) = dx(i + 1)
        dy(i + 2) = dx(i + 2)
        dy(i + 3) = dx(i + 3)
        dy(i + 4) = dx(i + 4)
        dy(i + 5) = dx(i + 5)
        dy(i + 6) = dx(i + 6)
   50 continue
      return
      end
      double precision function ddot(n,dx,incx,dy,incy)
c
c     forms the dot product of two vectors.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(1),dy(1),dtemp
      integer i,incx,incy,ix,iy,m,mp1,n
c
      ddot = 0.0d0
      dtemp = 0.0d0
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dtemp = dtemp + dx(ix)*dy(iy)
        ix = ix + incx
        iy = iy + incy
   10 continue
      ddot = dtemp
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dtemp = dtemp + dx(i)*dy(i)
   30 continue
      if( n .lt. 5 ) go to 60
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +
     *   dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
   50 continue
   60 ddot = dtemp
      return
      end
      subroutine dpofa(a,lda,n,info)
      integer lda,n,info
      double precision a(lda,1)
c
c     dpofa factors a double precision symmetric positive definite
c     matrix.
c
c     dpofa is usually called by dpoco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c     (time for dpoco) = (1 + 18/n)*(time for dpofa) .
c
c     on entry
c
c        a       double precision(lda, n)
c                the symmetric matrix to be factored.  only the
c                diagonal and upper triangle are used.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c     on return
c
c        a       an upper triangular matrix  r  so that  a = trans(r)*r
c                where  trans(r)  is the transpose.
c                the strict lower triangle is unaltered.
c                if  info .ne. 0 , the factorization is not complete.
c
c        info    integer
c                = 0  for normal return.
c                = k  signals an error condition.  the leading minor
c                     of order  k  is not positive definite.
c
c     linpack.  this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas ddot
c     fortran sqrt
c
c     internal variables
c
      double precision ddot,t
      double precision s
      integer j,jm1,k
c     begin block with ...exits to 40
c
c
         do 30 j = 1, n
            info = j
            s = 0.0d0
            jm1 = j - 1
            if (jm1 .lt. 1) go to 20
            do 10 k = 1, jm1
               t = a(k,j) - ddot(k-1,a(1,k),1,a(1,j),1)
               t = t/a(k,k)
               a(k,j) = t
               s = s + t*t
   10       continue
   20       continue
            s = a(j,j) - s
c     ......exit
            if (s .le. 0.0d0) go to 40
            a(j,j) = sqrt(s)
   30    continue
         info = 0
   40 continue
      return
      end
      subroutine  dscal(n,da,dx,incx)
c
c     scales a vector by a constant.
c     uses unrolled loops for increment equal to one.
c     jack dongarra, linpack, 3/11/78.
c     modified 3/93 to return if incx .le. 0.
c
      double precision da,dx(1)
      integer i,incx,m,mp1,n,nincx
c
      if( n.le.0 .or. incx.le.0 )return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      nincx = n*incx
      do 10 i = 1,nincx,incx
        dx(i) = da*dx(i)
   10 continue
      return
c
c        code for increment equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dx(i) = da*dx(i)
   30 continue
      if( n .lt. 5 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dx(i) = da*dx(i)
        dx(i + 1) = da*dx(i + 1)
        dx(i + 2) = da*dx(i + 2)
        dx(i + 3) = da*dx(i + 3)
        dx(i + 4) = da*dx(i + 4)
   50 continue
      return
      end
      subroutine dtrsl(t,ldt,n,b,job,info)
      integer ldt,n,job,info
      double precision t(ldt,1),b(1)
c
c
c     dtrsl solves systems of the form
c
c                   t * x = b
c     or
c                   trans(t) * x = b
c
c     where t is a triangular matrix of order n. here trans(t)
c     denotes the transpose of the matrix t.
c
c     on entry
c
c         t         double precision(ldt,n)
c                   t contains the matrix of the system. the zero
c                   elements of the matrix are not referenced, and
c                   the corresponding elements of the array can be
c                   used to store other information.
c
c         ldt       integer
c                   ldt is the leading dimension of the array t.
c
c         n         integer
c                   n is the order of the system.
c
c         b         double precision(n).
c                   b contains the right hand side of the system.
c
c         job       integer
c                   job specifies what kind of system is to be solved.
c                   if job is
c
c                        00   solve t*x=b, t lower triangular,
c                        01   solve t*x=b, t upper triangular,
c                        10   solve trans(t)*x=b, t lower triangular,
c                        11   solve trans(t)*x=b, t upper triangular.
c
c     on return
c
c         b         b contains the solution, if info .eq. 0.
c                   otherwise b is unaltered.
c
c         info      integer
c                   info contains zero if the system is nonsingular.
c                   otherwise info contains the index of
c                   the first zero diagonal element of t.
c
c     linpack. this version dated 08/14/78 .
c     g. w. stewart, university of maryland, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,ddot
c     fortran mod
c
c     internal variables
c
      double precision ddot,temp
      integer case,j,jj
c
c     begin block permitting ...exits to 150
c
c        check for zero diagonal elements.
c
         do 10 info = 1, n
c     ......exit
            if (t(info,info) .eq. 0.0d0) go to 150
   10    continue
         info = 0
c
c        determine the task and go to it.
c
         case = 1
         if (mod(job,10) .ne. 0) case = 2
         if (mod(job,100)/10 .ne. 0) case = case + 2
         go to (20,50,80,110), case
c
c        solve t*x=b for t lower triangular
c
   20    continue
            b(1) = b(1)/t(1,1)
            if (n .lt. 2) go to 40
            do 30 j = 2, n
               temp = -b(j-1)
               call daxpy(n-j+1,temp,t(j,j-1),1,b(j),1)
               b(j) = b(j)/t(j,j)
   30       continue
   40       continue
         go to 140
c
c        solve t*x=b for t upper triangular.
c
   50    continue
            b(n) = b(n)/t(n,n)
            if (n .lt. 2) go to 70
            do 60 jj = 2, n
               j = n - jj + 1
               temp = -b(j+1)
               call daxpy(j,temp,t(1,j+1),1,b(1),1)
               b(j) = b(j)/t(j,j)
   60       continue
   70       continue
         go to 140
c
c        solve trans(t)*x=b for t lower triangular.
c
   80    continue
            b(n) = b(n)/t(n,n)
            if (n .lt. 2) go to 100
            do 90 jj = 2, n
               j = n - jj + 1
               b(j) = b(j) - ddot(jj-1,t(j+1,j),1,b(j+1),1)
               b(j) = b(j)/t(j,j)
   90       continue
  100       continue
         go to 140
c
c        solve trans(t)*x=b for t upper triangular.
c
  110    continue
            b(1) = b(1)/t(1,1)
            if (n .lt. 2) go to 130
            do 120 j = 2, n
               b(j) = b(j) - ddot(j-1,t(1,j),1,b(1),1)
               b(j) = b(j)/t(j,j)
  120       continue
  130       continue
  140    continue
  150 continue
      return
      end
